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ABSTRACT 


Factorization is an approach to linear programming (LP) in which the algebraic 
elements of the LP tableau are organized in such a way that a large portion of the 
tableau may be represented implicitly and generated from the remaining explicit 
part. In dynamic row factorization, the row structure of the LP model instance 
influences the algebraic structure of the tableau, and the dimension of the algebraic 


elements may change as the solution progresses. 


We present three algorithms motivated by this approach, each resulting from 
a different LP model row structure: generalized upper bound (GUB) rows, pure net- 
work rows and generalized network rows. We describe implementations of all three 
algorithms, specifying data structures for tableau and basis inverse representations 


and detailing procedures for manipulation and update of these representations. 


Computational results are presented for a number of real-world models taken 
from a variety of applications and industries. From each model, one or more partic- 
ular instances are solved by each of our three implementations and by a commercial- 
quality mathematical programming system. The characteristics of the four solvers 
are compared and contrasted. Previous research on related algorithms by others 
suggests that these algorithms are properly viewed as specialized approaches, useful 
only on narrow classes of problems. Our computational results strongly refute this 
view, and instead suggest that each aj em is superior to the general simplex 


approach on a wide range of problem classes and structures. 
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I. INTRODUCTION 


A. INTRODUCTION 


А recurring theme in the development of simplex-based algorithms for linear 
programming has been the identification and exploitation of special problem struc- 
ture. Ideas as apparently disparate as the simplex method for bounded variables, 
primal and dual decomposition methods, pure and generalized network primal sim- 
plex algorithms and primal partitioning schemes may be unified to a degree by 
interpreting their development in this context. 

The factorization approach due to Graves and McBride [1976] provides a unify- 
ing framework from which we may reinterpret many existing algorithms and, through 
the application of the common principles embodied in the approach, develop new 
algorithms. This approach has as its central thesis the idea of recognizing, isolat- 
ing and exploiting special structures which may occur in a certain type of linear 
programming tableau. An example of such special structure is generalized network 
rows. Generalized network rows are naturally specified as a row structure in which 
each column has at most two nonzero elements within the rows. 

This paper adopts the design principles and algorithmic structure suggested 
by Graves and McBride [1976] to develop algorithms that, while general in nature in 
the sense that each may be used to solve any linear programming problem instance, 
are strongly tailored to exploit a particular row structure. We call this approach 
“dynamic row factorization”; “row factorization” because we exploit the structure 
in the basic tableau which is induced by the row structure of the LP model instance, 


and “dynamic” because the dimension of the structure may vary (or even fail to be 


present) as the solution progresses. In our setting, we require the row structure of 
the model instance to be specified prior to solving, and to remain fixed throughout 
the solution process. Ап extension of this approach is to allow the row structure of 
the model instance to vary as the problem is solved. While we do not develop an 
implementation of the second approach, we show that it is a conceptually simple 
extension of our work. 

Each algorithm is developed by factoring the constraints of the LP model 
into two classes: those that have special structure (factored) and those that do 
not (explicit). This factoring of constraints induces a factored structure in the LP 
tableaus which may be exploited computationally. We treat each of three structures 
of factored constraints in this work: generalized upper bounds (GUB), pure networks 
and generalized networks. 

We implement each of the three algorithms by integrating it within the struc- 
ture of a state-of-the-art mathematical programming system: the X-System of 
Brown and Graves [1975]. We do so both to demonstrate the feasibility of imple- 
menting such methods and to determine whether or not such methods have practical 
value for the practitioner interested in solving real-world problems. Commercial im- 
plementations of similar GUB algorithms have generally failed to inspire enthusiasm 
and are apparently falling into disuse [Kennington [1978]]. While the first commer- 
cial implementation of a related “pure network with side constraints” algorithm is 
(to our knowledge) just now becoming available, reports of research implementations 
have generally supported the view that such algorithms have limited application. 
The reports of implementations of similar “generalized network with side constraint” 
algorithms have offered less promise than their pure network counterparts. 

After we review the related literature, we provide a detailed presentation of the 


underlying tableau-based algorithm which forms the foundation of all our subsequent 


work. We do so because the algorithm is not widely known and may be unfamiliar 
to the reader. We then review the general factorization approach of Graves and 
McBride [1976]. We provide a design template for our developmental approach, and 
then present the dynamic row factorization algorithms for GUB, pure network and 
generalized network row structures. Computational results are presented, and we 


then summarize our conclusions and suggest avenues for further research. 


B. SURVEY OF RELATED LITERATURE 

While the terms “partitioning” and “factorization” are frequently used inter- 
changeably in the literature, we observe a distinction between the two approaches. 
We consider partitioning methods to be those that are based on special structure 
in the original problem instance. This structure in the problem instance need not 
induce special structure into the LP tableau, and in fact the method need not be 
tableau-based. This is in contrast to our view of factorization, in which the algo- 
rithm is based on special structure which occurs in the simplex basis and thus in the 
basic tableau. We thus classify Dantzig- Wolfe [1960] and Benders [1962] Decompo- 
sitions and Rosen’s [1964] primal partitioning method as examples of partitioning 
methods. 

Perhaps the first example of what we consider factorization 1s the treatment of 
simple upper bounds by Dantzig [1954] and, independently, by Charnes and Lemke 
[1954]. They observe that it is more efficient to enforce the “structural” simple 
upper bound constraints with logical tests within the algorithm rather than treat 
them explicitly along with other constraints. While not originally presented in the 
context of a formal tableau factorization. the approach is easily viewed as such and 


is consistent with the general approach. 


The mutual primal-dual method of Graves [1965] focuses attention on the 
special role of nonnegativity constraints in linear programming. A clear distinction is 
drawn between the computational convenience of treating nonnegativity constraints 
implicitly rather than explicitly and the unambiguous mathematical equivalence 
of all problem constraints, structural or nonnegativity. Emphasizing the special 
importance of inequality constraints, the approach yields an elegant theory (see 
Graves [1987]) and, as we will see, efficient implementations. We view this algorithm 
as the first formal example of factorization. 

Dantzig and Van Slyke [1967] extend the factorization approach applied earlier 
to simple upper bounds in a more structured manner in their treatment of general- 
ized upper bounds (GUB). In a problem with p GUB constraints and m structural 
constraints, their approach requires a working basis inverse of dimension (m+ 1), a 
considerable savings when p is large. 

Hartman and Lasdon [1972] specialized this approach to the multicommodity 
capacitated transshipment problem. In this case, the structure of the basic pure 
network columns introduces additional structure into the working basis, allowing 
further simplifications in basis representation and update techniques. Helgason and 
Kennington [1977] develop techniques for representing the working basis inverse in 
product form and provide graphic interpretation of the graph updates. Kennington 
[1977] reports an implementation of the algorithm. 

McBride [1972] and Graves and McBride [1976] formalize and generalize the 
factorization approach. They view it is as a unifying framework for tableau-based 
simplex specializations and illustrate this by developing a variation of the GUB 
algorithm of Dantzig and Van Slyke [1967] and a GUB algorithm for doubly coupled 
linear programs of Hartman and Lasdon [1970]. They present a new algorithm 


for the set partitioning linear programming problem and an equality-constrained 


form of the pure network with side constraints model. McBride [1972] reports ап 
implementation of a GUB factorization. 

Schrage [1975] extends the approach of simple and generalized upper bounds 
by considering variable upper bounds (VUB), which are constraints of the form 
r; S Tk, Where z, is said to be the variable upper bound of z;. His algorithm allows 
the implicit representation of the VUB constraints by expressing VUB variables in 
terms of other basic columns. This permits the basis representation to be treated 
in two parts, one part a large matrix which changes infrequently and thus needs 
to be updated only occasionally, and the second part a small working basis which 
requires regular attention. Thus, computation and storage savings may be realized. 
Schrage [1978] extends these ideas to what he calls generalized variable upper bounds 
(GVUB) constraints, which arise frequently in models involving fixed charges. 

Klingman and Russell [1975] sketch a factorization method for solving trans- 
portation problems with side constraints. They suggest techniques for performing 
simplex iterations and updating the problem representation. Chen and Saigal [1977] 
present a similar approach for solving capacitated network flow problems with ad- 
ditional linear constraints. Both the above presentations consider a graph-theoretic 
view of the basis update mechanism and allow the basis representation to be treated 
in two parts, a part which corresponds to a rooted spanning tree defined on the un- 
derlying graph, and a general working basis inverse. Glover, Karney, Klingman and 
Russell [1978] report an implementation of the Klingman and Russell design, but 
one which (curiously) only accommodates a single side constraint. McBride [1989] 
reports an implementation which requires the pure network rows to be equalities 


and allows more than one side constraint. 


The problem of generalized network problems with side constraints is ad- 
dressed by Hultz and Klingman [1976]. They present details for the simplex price- 
out, column generation and basis update. Hultz and Klingman [1978] report an 
implementation that (curiously) solves the “singularly constrained” generalized net- 
work problem. McBride [1989] reports an implementation that is not restricted to - 
a single side constraint. 

The factorization approach has been extended by the consideration of em- 
bedded structures. Glover and Klingman [1981] consider a linear program which 
contains embedded pure network structure, 1.e., the pure network structure appears 
in only a subset of the rows and columns of the technological coefficient matrix of 
the problem. Their approach yields an algorithm similar in spirit to the algorithms 
for the pure network with side constraint model, but the presence of the “side vari- 
ables" significantly complicates the basis representation and update. They report 
an implementation of the algorithm but (curiously) restrict the problem suite to 
problems having no complicating variables. 

McBride [1985] treats the problem of linear programs with embedded general- 
ized network structure. He presents methods for pricing, column generation, basis 
representation update and data structures. A successful implementation is reported 
which is approximately five times faster than MINOS [1977] for the models tested. 

Interest in developing algorithms to solve problems with special substructures 
has been accompanied by work to identify such substructures in problem instances. 
Greenberg and Rarick [1974] and Brown and Thomen [1980] develop algorithms 
to identify GUB sets. Brown and Wright [1984] develop algorithms for identifying 
pure network constraint substructures. Brown, McBride and Wood [1985] present 
a method for locating generalized network structures, both embedded and row-only 


structures. 
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Todd [1983] examines factorization from a geometric standpoint and constructs 
a geometric interpretation which is in large measure equivalent to the algebraic 


development of Graves and McBride [1976]. 


II. MUTUAL PRIMAL-DUAL METHOD 


A. INTRODUCTION 

This chapter presents the mutual primal-dual linear programming method in- 
troduced by Graves [1965] which provides the algorithmic framework and notational 
conventions for the research which follows. 

We begin with a complete algebraic development of the primal algorithm fol- 
lowed by a less detailed symmetric discussion of the corresponding dual algorithm. 
We conclude with a unified treatment of the two which establishes the theoretical 
importance of the algorithm and justifies its use as the foundation of our specializa- 
tions. 

The following presentation scrupulously restates the Graves [1965] algorithm, 
incorporates later discussion by McBride [1972], and accommodates large-scale im- 
plementations by Brown and Graves [1975]. The view presented here is not available 


in standard reference texts, and is included in the interest of completeness. 


B. PRIMAL PROBLEM STATEMENT (PLP) 


The traditional statement of the linear programming (LP) problem is: 


(LP) min: wy 
у 
s du S Т Е 
еу 2 0 ‚ Ј= 1,...„п, 
where y is an n-vector of decision variables, w is an n-vector of cost coefficients, 
each a; is an n-vector of technological transformation coefficients, each r; is a scalar 


right-hand side coefficient, and e; is the j'^ unit vector. While this statement of 


the problem is clear and unambiguous, there are reasons for preferring an alter- 
native. The insistence upon drawing a formal distinction between the "structural" 
constraints a;y € r; and the “nonnegativity” constraints e;y > 0 obscures the math- 
ematical structure of the problem by suggesting that the two types of constraints 
are inherently different. Certainly the exploitation of the special structure of the 
e;y > 0 constraints leads to computational efficiencies in the implementation of the 
algorithm. However, in the theoretical development of the algorithm, we prefer to 
treat them simply as general inequality constraints. 

In order to achieve a consistent form, we rewrite the nonnegativity constraints 
as —e;y € 0 and group them with the structural constraints. The problem statement 
then becomes: 

(РЛЕР) min : шу 


See aly E T; „0 тя 


where wy is called the extremal function. The constraints of (PLP) define a set of 


feasible solutions which we shall call the feasible set, F: 


Ее {ує R"|ay&ri, Ее ооа 


Since F is the intersection of a finite number of closed half-spaces, F' is itself a closed 
set. If F is nonempty and bounded, then an optimal solution to (PLP) will occur 
at an extreme point of F. 

A point y? € F is said to be a feasible point or feasible solution. If, for 
constraint 7, a;y° <r, . constraint i is satisfied at y? , and the quantity r; — a;y? is 
the slack in constraint i at y? . If, on the other hand, for constraint 1, ay? » ri , 
constraint i is violated at the point y? , the magnitude of the violation being a;y? — r, 


(the negative of the slack). 


A point y? € R” is defined to be a basic solution of (PLP) if there exists an 
independent subset {a;,,@;,,.-.,@i,} Of {@1,42,.--,@m4n} such that aj,y? — ri, for 
j =1,...n. Such an independent subset {a;, , aiz,- .., Qin} of {a1,@2,...,4m4n} is a 
basis for R” , and y° is the (basic) solution to the system a; y? = т;,,ј = 1,...,п. 

Each basic solution of (PLP) corresponds to an eztreme pointof F and for each 
extreme point of F there is at least one corresponding basic solution of (PLP). Since 
there are at most ( rAr ) ways of choosing an independent subset of n vectors from 
{a,,42,...,@m4n}, the number of basic solutions of (PLP) and thus the number of 
extreme points of F is finite. Hence, (PLP) can be solved by searching among its 
basic solutions. The algorithm to be developed here will implicitly enumerate the 


set of basic solutions of (PLP) and terminate in one of three states: 


IE = ( 


(no feasible solution exists) ; 


2. the extremum is unbounded 


(for every real number a there is a point y? € F with wy? « o ); or 


J. there exists at least one optimal solution 


(a point y" € F with vy" wy Vy € F ). 


We will first consider the problem of finding a basic feasible solution to (PLP). 
Having achieved this, we will then consider the task of finding a basic feasible solution 


which is also optimal. 


C. OBTAINING FEASIBILITY 
Since we have included the nonnegativity constraints —e;y <0, 7 =1,...,n 
in our structural constraints ayy <т7,,і-1,....т--п, and since the origin is the 


unique solution to the independent system: 
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—е;у = 0, J= binas 


the origin is a basic solution that is always immediately available for (PLP). 

Let y? be any basic solution to (PLP). By definition, there are at least n linearly 
independent constraints which are exactly binding at y? . Among the m remaining 
constraints, typically some will be satisfied at y? and others will be violated. Our 
strategy will be to focus on one violated constraint at a time, which we will call 
the target constraint. Moving from one basic solution to another, we will attempt 
to reduce the violation of the target constraint until it becomes satisfied. We will 
restrict our choices of basic solutions in such a way that all constraints that are 
already satisfied at y? remain satisfied at each subsequent basic solution. Once the 
target constraint becomes satisfied, we then select some other violated constraint 
as the new target constraint, and repeat the process. We proceed until either all 
constraints are satisfied and we have obtained a basic feasible solution, or we find 
a violated constraint which cannot be satisfied, in which case we conclude that no 
basic feasible solution exists. 

To formalize these ideas, define S(y°) to be the set of indices of all constraints 


satisfied at a basic solution y? : 


500%) = {1 <і<т+п [ау < т). 


Of course, | $(y9) | > п. 
Suppose constraint k is violated at the basic solution y°. Then axy? > ry, and 
k d S(y9). A necessary condition for the existence of a basic feasible solution is that 


there exists a basic solution y with: 


TH 


5(у) > 5(у?) (2.1) 


and 


aky < aky”. (2.2) 


If there exists no such basic solution y satisfying (Eq. 2.1) and (Eq. 2.2 ) we conclude 
that F — 0. Thus, we may restrict our attention to basic solutions satisfying (2.1) 
and (2.2). 

[её {а;,,а;,,...,а;,} Бе а basis for E" at y?. For notational convenience we 
will partition the constraints into two sets, those that are basic at y? and those that 


are nonbasic at y: 


bi ai, Л ri 
Е | аа (2.3) 
b, Gig f. d 
di а, gi ШТ 
ee Mm ‚ в#ө=| fal (2.4) 
d ТАРА 9т е 


Define B and D to be the set of indices of the rows of B (basic constraints) 
and D (nonbasic constraints) respectively. Using this notation, the current basic 
solution y? may be expressed as By? = f, and since the rows of B are by definition 
linearly independent, B^! exists and y? = B“! f. 


Let y be any other point of R”. Then y can be expressed as: 


у= у + (у = у). 


We have chosen this representation since at y? the basic constraints are satisfied 
exactly and thus it is the direction vector (y — y?) of any proposed move from y? 
to y that determines whether or not a given basic constraint will remain satisfied 
at y. ( Of course, it is the magnitude of (y — y?) that is important in determining 
whether the nonbasic constraints rernain satisfied as well.) 

Now, for a given arbitrary basis (pl,p?,..., p^) of R” and any point y, there 


must exist scalars Àj, ..., A, such that: 


у=у°+(у—у°) = у°+ У Asp’ =y + PA 


2-1 
Since the choice of P is arbitrary, let us choose P — —B^!, which we call the 


conjugate row basis. Then: 


Ву = В(у? + РА) = By’ — BA) = By! - AS f — A, 


and thus y satisfies the basic constraints By € f if and only if A > 0. 
Thus, the set of points in R” which satisfy the basic constraints b;y € fj for 


1 Є B can be characterized as: 


ен ya 4D awd, 2078) 


JEB 
“ith P = - B7! . 


It then follows that everv point y satisfying 


E ME (2.5) 
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can be expressed as 


у= у + РЛ, 


for some À » 0 апа Р = – В-!. 
To address condition (2.2), assume that we have selected constraint k as our 
target constraint (and thus k is violated at y? , i.e., dy? > gk). To reduce the 


violation of the target constraint, we seek a point y such that: 


dky < dky’, (2.6) 


holds. Since such a point must also satisfy (2.1), it must have a representation of 


the form y = y? + PA for some à > 0. Replacing y by y° + PX , we have: 


dily? + PA) < ау, 


which holds if and only if d, PA < 0. Since à > 0, it follows that a necessary 
condition for the existence of a point y satisfying (2.5) and (2.6) is that at least 
one element of the vector d,P be negative, i.e., d,p’ < 0 for some j, 1 <j <n. If 
d, P > 0, we conclude that F = @. 


Suppose that dyp' < 0. Then (2.6) holds at all points of the form: 


1 = yo + Ap’, Ar 06 


The generic point y? -- Ajp! lies along an edge of the set of all points satisfying the 
constraints which are basic at y? . If there is more than one l for which dp! < 0, 
any one can be chosen. For the purpose of exposition, we will designate the first 


such index encountered as “l”. 
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The violation in constraint k decreases linearly as A; increases. Constraint k 


becomes satisfied when: 


дь(у° + Мр) = gi, 


or when à, = t, with 


A (gk — dky?) 
t — I—————ÀÀ А 
d, p! (2.7) 


А geometric interpretation is that y = y° + tp’ is the point at which the ray 
{у Riy =y t Ap, À 0} pierces the hyperplane {y E€ R” | dky = gx}. 
Define àf to be our ultimate choice for 4; . Choosing t as the value for Aj 
satisfies the target constraint k, but we are required by (2.5) to continue to satisfy 
all nonbasic constraints which are already satisfied at y? (if any). The choice of t 


may cause violations of such constraints. Thus, À? must also satisfy the condition: 


d;(y° + Xp')<9; Y j EDAS’). 


Writing this as 


Дар <4;- d;y, 


we find that we must choose A, € s, where 


(om d;y^) 
2-5 а; 0. 2.8 


If this set is empty, define s = œ . Of course, it is possible for s to be equal to zero. 


A 
ET 


Thus. by choosing 


А 2 min [146 |; (2.9) 


where t and s are as defined in (2.7) and (2.8) respectively, we obtain the largest 
reduction in the violation of the target constraint consistent with the feasibility 
restriction (2.5) with respect to the chosen direction p! . 

The selection A? according to (2.9) leads to a new basic solution y! — y? -- Ap p/,. 
If A7 » 0 at each step, the transitivity of (2.5) and (2.6) guarantees that no basic 
solution will be repeated before any given target constraint is satisfied, or until 
the conclusion is reached that F = 0. At any given step, either the feasible set is 
shown to be empty, a basic solution is found which satisfies the target constraint 
or a positive reduction in the violation of the target constraint is achieved as the 
result of moving to another basic solution. Since the set of basic solutions is finite, 
the third alternative may be repeated a finite number of times before one of the 
first two alternatives occurs. If this constraint becomes satisfied, (2.5) guarantees 
that every subsequent solution will also satisfy the constraint. We may then select 
a new violated nonbasic constraint as the target constraint. Since there are a finite 
number of constraints, we will either discover a basic feasible solution or determine 
in a finite number of steps that none exists. 

If Ay = 0 at any step of the algorithm, we say that a block has occurred. 


Blocks will be discussed later in detail. 


D. FROM FEASIBILITY TO OPTIMALITY 
Once feasibility is achieved, say at y?, the process of proceeding to optimality 
can be thought of as minimizing over F the violation of the constraint wy € wy?— a, 


where a 1s a sufficiently large positive constant. 
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Let y? be a basic feasible solution with By? — f and P — —B-! . Since every 
point y E€ R” which satisfies By € f may be written as y — y? -- PA for some \ > 0, 


the value of the extremal function wy at such a point y is: 


шу = w(y? + PA) = wy? + wPA = wy? + 9 лушу. 


j=l 
Thus, a necessary condition for the existence of a point y € F such that wy € wy? 


15: 


шр? < 0 for some ПБ п. 


If wp? 2 0, then y? must be a feasible, optimal solution to (PLP). 

Suppose wp! < 0. Since all constraints are satisfied at y°, the greatest reduc- 
tion in the value of the extremal function in the direction p' is achieved when A7 is 
chosen to be s, where s is as defined in (2.8). If s = œo , (PLP) has an unbounded 
extremum. If 0 « A7 « oo , a positive reduction in the value of the extremal function 
can be achieved by moving from the basic feasible solution y? to the basic feasible 
solution y = y? + A;p! . If A? = 0, a block is encountered. 

Once a feasible solution has been obtained, (2.5) and (2.6) become equivalent 


to: 


y€F , (2.10) 


wy < wy? . (2988) 


The transitivity of (2.10) and (2.11) ensures that no basic feasible solution will 


be repeated as long as A7 » 0 at each step. At each iteration, either a basic feasible 


ІН 


solution is found to be optimal, the solution is found to be unbounded or a positive 
reduction in the value of the extremal function 1s achieved by moving to another 
basic feasible solution. Since the set of basic feasible solutions is finite and no basic 
feasible solution is repeated, the third alternative can only occur a finite number of 
times. Thus, if A7 » 0 at each step and if a finite optimal solution exists, it will be 


found in a finite number of steps. 


E. BLOCKS 


Suppose that y? is a basic solution satisfying 


where, as before, B is of dimension n by n and nonsingular and p! is the j'^ column 
of P 2 — B^!. Then, y? is the unique point in R” lying in the intersection of the n 


hyperplanes: 


осыда а 


Suppose that p’ is a direction that either leads to a reduction in the violation 
of some target constraint or, if y? is feasible, a reduction in the value of the extremal 


function. Since the nonbasic constraints are of the form: 


a block occurs when. for at least one nonbasic constraint k, 


diy" = gi. 


and 


Фр! >0. (2.12) 


Thus, any movement away from y? in the direction p! will result in violation of the 
nonbasic constraint d,y9 < gk . 

If dky? = gk , then y° is in a sense *over-determined", and y? is said to be 
a "degenerate" solution. Geometrically, y? lies in the intersection of at least n + 1 
hyperplanes. Algebraically, there is more than one basis that can be formed from 


the row vectors {a}, @2,...,@m+4n} for which: 


D m 
Qi, y se 3 J т, 


holds, and there is more than one basic solution which corresponds to the extreme 
point y? . A block is therefore encountered when an over-determined solution sat- 
isfies (2.12) in an “improving” direction (e.g., one that leads either to a reduction 
in infeasibility or. if feasible, to a reduction in the value of the extremal function). 


We will develop a method for dealing with blocks later. 


F. THE PRIMAL ALGORITHM AND A SUPPORTING BASIC TAB- 
LEAU 


The primal algorithm proceeds as follows: 


1. Identify an initial basic solution. Notice that the origin satisfies 
le y^ ES 0, 


and thus mav alwavs serve as the initial solution. 
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2. If the current solution is infeasible, select a violated primal constraint index k 


(i.e., an index k for which g, — d,y? « 0). This requires the quantity: 


g-Dy . 


3. Within the target row, select an element of the proper sign (negative, according 
to our convention) whose index, l, specifies a “transformation column”. If the 


current solution is infeasible, this requires the quantity: 
Or 


with the transformation column satisfying dąp' < 0. If the current solution is 


feasible, this requires the quantity: 
wP 


with the transformation column satisfying wp! « 0. If no such element exists, 
then the problem is infeasible (if the current solution is infeasible and d, P > 0) 
or optimal (if the current solution is feasible and wP > 0). The task of 
selecting such an index l is commonly referred to as “pricing” or a “pricing 


strategy”. 


4. Compute t as in (Equation 2.7). If the current solution is feasible, assign 


= Ес 


9. Compute s as in (Equation 2.8). This computation is commonly called a “ratio 


test’, or a “minimum ratio test”. This computation requires the quantities: 


g— Dy? and Dy. 


6. Compute A; as in (Equation 2.9). If A7 = ос, the problem is unbounded. 
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7. Update the current solution according to the computation: 
у= у AP. 


8. Update the assignments of constraint indices to B and D, and update P. 


9. Go to Step 2. 


The only quantities needed at each step of the algorithm are the matrix DP, 
the column vector g — Dy? and the row vector wP. For notational convenience, we 
can define D to be an (m 4- 1) by n matrix whose bottom row dm41 is w, and 4 %о 


be a (m 4- 1) dimension column vector whose bottom entry, 4541, is zero. Then 


Gmt+1 — dingy” -0- шу° = —wy^, 


which is the negative of the extremal function value at the point y . We can now 
conveniently display all the relevant information by forming the (m -- 1) by (m 4- 1) 


matrix: 


[DP] ġ- Dy], (2.13) 


which we will call the basic tableau. 

By computing the basic tableau in partitioned matrix form, we may isolate the 
important algebraic components required by this method. To do so, let us assume 
that at the current basic solution the basis consists of h structural constraints and 


(n — h) nonnegativity constraints. Then (2.3) can be written as: 


21 


bi а, fi Ti 
bz au h m 
b, а; Sh Ti 
B= = B and f — == a 
брт Cina f Ман 0 
bn+2 cc Аһ? 0 
Б, — Cjn Ja 0 
In partitioned matrix form, 
k n—h 
и >-- 
h 
B = Ап А12 } (2.14) 
0 -І)/ )n-h , 
and thus 
h n—h 
~ ~ 
р 8 Ed —Ai A ) h 
0 1 }n—h 
Similarly, (2.4) can be written as: 
di € д1 0 
d — 3, 92 0 
" M 
Р = P ES Ejn and Qum Sh = 
А+1 а; вт 9һ+1 Tini 
dj 42 Ging? 9h42 Ving? 
dm Com Gm Dim 


and in partitioned matrix form: 


и ,-- 
= h 
же (2.15) 
Ал А, үрт == h г 
Then 
À n—h 
=“ a 
Ат Aj Ai } h 


DP = -DB^ = 
—АА у Ao — АлАТАр/ |т-һҺ, 


and we shall call DP the principal part of the tableau. By partitioning w = 
(wi, w2), 91 — (9,92)! and y? — (y9, y9), the basic tableau (2.13) may be written 


in partitioned matrix form as: 


к = A. 

Ат Ay Ai у? )h 
—AnÁi А)- АлАПА) g2— Any? — Ar2ys | )m- h (2.16) 
-wAn w- wAn A? —wy? Ы. 


Note that y? is displayed explicitly in the tableau (2.16). Also, y9 — 0 since the 
corresponding nonnegativity constraints are basic and thus binding. 

As the algorithm proceeds from one basic solution to another, (2.16) can be 
updated using a slight variant of the Gauss-Jordan transformation (which we will 
call a pivot). The transformation is as follows: 

Let |D P|g — бу = [v,;] be the basic tableau at any basic solution y? . If 


vs, = d,p' # 0 then the basic tableau at the basic solution: 


1__ 0 Us.n41 lx 0 (gs — dsy°) l 
U 0 Ел = у +(e) үн 
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written as 


|D'P' | 4! Ке D'y' = lvii 


is calculated as: 


general elements (1 s,j x 1l): 


Us. 504.1 
/ — No Í — R57 n 
©» 1 


pivotal column (2955,7): 


ке — Uil 
йр? , 
9,1 


pivotal row (nes s EE 


/ == Us,j 
Cr es 
J 
951 


pivotal element ғ: = $, у = і 


d — 
9,1 a 





v 

After applying the transformation to (2.16), row and column interchanges may 
be necessary to restore the explicit structure of B in (2.14) and D in (2.15). A proof 
of this may be found in Graves [1987], who also extends the basic algorithm to 


include free variables. equality constraints and bounded variables. 


G. DUAL LINEAR PROBLEM (DLP) 


Corresponding to (PLP) is the following linear program called its dual: 


(DIET max: тт 
st: TAK W 
si- o 


The relationship between (DLP) and (PLP) is as follows: 
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1. if (PLP) is infeasible, then (DLP) is either infeasible or unbounded. 
2. if (PLP) is unbounded, then (DLP) is infeasible. 


3. if (PLP) possesses a finite optimal solution y* , then (DLP) also possesses a 


finite optimal solution z* and wy" — z*r. 


See Graves [1987] for a proof of this. 

Solving (DLP) directly has important advantages in certain problems, and we 
will show that its solution plays an important role in dealing with blocks in (PLP). 
We will thus develop a dual algorithm in a manner similar to that for (PLP), sparing 
some of the symmetrical detail where possible. 

Гей т, Бе а basic solution for (DLP). Then there exists an independent subset 
IN (5 ор (0142... 2777! such that тоа — w^, i = 1,...,m, and 
this linearly independent subset will be called the basis for R" at z, . 

Partitioning the dual constraints into those which are basic at z, and those 


which are nonbasic at z, yields: 


T = (Аа... Кезі! = (ал ал)... а:”) 


и = (u',u?,...,u™) = (w^, w?,...,w!7), (21) 
апд 
О РЕ (ао ат) 
v= (utes... et) = (wit, wit? "). (2.18) 


Define 7 and X to be the set of indices of the columns of T (basic dual 
constraints) and A (nonbasic dual constraints), respectively. The current basic 


solution ze may then be expressed as 


to = uT}, 


and any other point z in R™ may be expressed as 


стат ы Ta 


For a given arbitrary basis {q),92,---;@m} of R™ , any such point z may be repre- 


sented as: 


r= х, + У" = To + YQ. 


t=] 


The current basic constraints are ra^ — w^,i — 1,...,m and for these con- 
straints to remain satisfied at the generic point zt = zo + YQ, we must have 
[ra^ € w^,i 2 1,...,m]. Choosing Q — T^! as a convenient basis for R™ at z,, we 
then have: 


tT = (zat YQ)T = trof ОСЕТЕ 


and thus we must have v» < 0. 


Now consider a violated dual nonbasic constraint t, where a violation means 


rok cw. 


To reduce the violation, we seek a point z at which 


tk e rk 
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holds. Such a point must also satisfy the basic constraints, and thus the condition 


for a reduction in the violation of the dual target constraint is: 


ie (few kh ГІЗІ ТІ 


which implies that 


ХОК" < 0. 


Since y% < 0, a necessary condition for constraint t to be satisfied is that there exists 
ап 1 with q;k' > 0. If none of the elements of q;k' are positive, we conclude that the 
dual constraints are inconsistent. 


The dual algorithm proceeds along the edges 


- А. 
T= Te 01. 


from basic solution to basic solution. We define 


S(z.) = {1 Sj S m+n | zai € ui], 


and insist that 


Sarl О (тж 


To satisfv this condition, we consider the effect of moving along the edge of a 


general nonbasic dual constraint k? which is currently satisfied at ro. We must have 


tO 
--1 


т? = (1, +) = т, + < 07. 


Since the constraint k? is satisfied at z, , we have 


toki < v > v — zok! > 0, 


and since v! « 0, if q;k? « 0 then as v)! decreases in value we approach the boundary 
of the constraint. Ё 02 — z,k? = 0 then q;k? < 0 causes an immediate violation if 
V* « 0. Thus, a block has occurred. 


We have shown that our choice for 2+, denoted y', should be 


p = max {b,c}, (2.19) 
where 
A qt t t 
b= (vu = rok yoke (2.20) 
and 
c S max { (0? — zk?) /q;k? |v — 2k? > 0, gik? < 0) (2:219) 
зек о 1 о ЖЕЕ Т % ч % 


] 
Once dual feasibility has been achieved, we wish to maximize: 


vr ae a aa 


Since uv < 0, a necessary condition to achieve an increase in the value of the extremal 
function is that there exists an ? such that qır « 0. Hence, when q;r 2 0, we conclude 
that т„ 15 dual optimal. 


The dual algorithm then proceeds as follows: 


2 
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1. Identify an initial basic solution. Notice that the origin satisfies 
29.1] = 0, 
апа thus may always serve as the initial solution. 


2. If the current solution is infeasible, select a violated dual constraint index t 


(1.e., an index t for which v' — zok' « 0). This requires the quantity: 
Uem тол. 


Column t is then referred to as the "target column". If the current solution is 


feasible, the right-hand side column is designated the “target column”. 


3. Within the target column, select an element of the proper sign (positive, ac- 
cording to our convention) whose index, ?, specifies a “transformation row”. 


If the current solution is infeasible, this requires the quantity: 


Qk , 


with the transformation row satisfying g;k' > 0. If the current solution is 


feasible. this requires the quantity: 


От › 


with the transformation row satisfying g;r > 0. If no such element exists, then 
the problem is infeasible (if the current solution is infeasible and q;k* < 0) or 
optimal (if the current solution is feasible and q;r € 0). The task of selecting 


such an index 7 is commonly referred to as “pricing”, or a “pricing strategy’. 


4. Compute k as in (Equation 2.20). If the current solution is feasible, assign 


b 2 —oc. 


5. Compute c as in (Equation 2.21). This computation is commonly called a 


*ratio test". This computation requires the quantities: 
v— r9K апа ОК. 
6. Compute v? as in (Equation 2.19). If i = 00, the problem is unbounded. 
1. Update the current solution according to the computation: 
їі = To + V q;. 
8. Update the assignments of constraint indices to 7 and K, and update Q. 


9. Go to Step 2. 


This development shows that at each step, the quantities QK, (v — z;, K) and 
Qr are needed to proceed with the dual algorithm. 

To develop a matrix partitioned form of the dual tableau, we proceed as before. 
Assuming the dual basis consists of À structural constraints and m — À nonnegativity 


constraints, we have: 


Т z(U,.,...,0) (ato NIE 


а (о О ЕУ 
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Emo С 1 
Ue es ооа 
The nonbasic constraints are then: 
a КИТТЕ | | | 
Е ЕО 2 


vm (vl... v^) 2 (0.0,.... 0 u7^*t1 u^, , u^). 


“ 


БЮ -- (0197) 2 (0,w?) . 


The matrix partitioned form of the basis is then 


À m—h 
езь« SN 
Ап 0 |А , 


and with the choice of Q = T~! 


h m—h 
ж- ~ 
ae Ai 0 | J^, 
—AnAy I mm 
with the remaining constraints forming 
‚ __ І А2 
h= | UO | 
The principal part of the dual tableau is: 
Р = Ay Of] 7 Ау; 
ae | VEN | | Oma | 
= Ат АТ А) 
—AnÁán А – АА А2 |, 


which we find to be exactly the principal part of the primal tableau, so 


I) P ru 
The quantities (v — z, A) are 


DEUM — t—(uQ)h 
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— v—uDP 


Asi ASA 
_ 2 11 11 4312 N 
= (0, w ) | 0 1 | uDP 
— vP-—uDP 
= (v—uD)P 


= (6.42) — (w!, 0) | " a ІІ 
= {(0,w?) —(—w’,0)}P 
cM 

The quantity Qr is: 


е = еп) 
- el(?)*(5.)] 
- e 21 (2)] 
= Q{94+ Kf} 
= Qg+QKf 


1 ТЕ 
= iig Hs) nor 


- (* рм 
= gra 


= g—Dy’. 


We thus find that the quantities required for the dual tableau (QK,v—z,K and 
Qr) are exactly the quantities required for the primal tableau (DP, wP and g— Dy). 
Therefore, the primal and dual algorithms use the same tableau. Reviewing the 


summaries of the primal and dual algorithms, we see the strong symmetry between 
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the two. Operating on a single tableau, the primal algorithm identifies a target 
row, selects a transformation column, performs a ratio test, and performs a pivot 
update, while the dual algorithm identifies a target column, selects a transformation 
row, performs a ratio test and performs a pivot update. This remarkable symmetry 
allows us not only to perform either algorithm on a single tableau, but to switch 
between one algorithm and the other as often as we choose. This fact has important 
implications for our implementation and it enables us to deal with blocking in a 


systematic and consistent way. 


H. RESOLUTION OF BLOCKING 


The resolution of blocking in either the primal or dual algorithm will be accom- 
plished by shifting to the alternate algorithm when blocking occurs. The alternate 
algorithm is applied to a subproblem of the original and at worst we are lead to a 
contracting sequence of problems to which we alternately apply the primal and dual 
algorithm. A strict contraction is assured, and thus in at most a finite number of 
steps resolution is achieved. 

We will demonstrate this procedure by assuming that we have started with 
the primal algorithm. When a block occurs, assume we have rearranged the rows 
so that the zeros in the right-hand side column occur contiguously. 

Let kj 2 n4 land k, — t where t is the index of the target row (if the current 
solution is primal feasible, kọ = m 4- 1, the extremal function row). Let k4 be the 
pivot column which caused the block. The situation is shown in Figure 2.1. 

Define Subproblem (1) as consisting of all the columns of the original problem, 
but only those rows with zeros in the right-hand side column, and as shown in Figure 
2.2. “extremal function” row kz. the target row of the original primal problem. Row 


ks is the row containing the blocking element. Now apply the dual algorithm to 
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Figure 2.2: Subproblem (1) 


Subproblem (1). Because the right-hand side column of any potential pivot row is 
zero, the right-hand side column of the original primal problem is invariant to any 
pivot in Subproblem (1). Let v;; denote the element in row i and column j of the 
tableau. We proceed with the dual algorithm on Subproblem (1) until one of the 


following situations occur (which may be after zero or more dual pivots): 


l. Ukk 2 0. A pivot in column Аз по longer reduces the violation of primal 
constraint k? , and thus the primal block has been eliminated. We thus return 


to the original primal problem and seek a new column in which to pivot. 


2. Uk, .k < 0. A pivot in column з по longer threatens to violate primal con- 
straint k4. We thus compute new ratios using columns k, and ky of the original 
primal tableau to determine if we can make a gain on the target constraint Ё, 


using column kz. 
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k2 


Ка 





k2 


Figure 2.4: Subproblem (2) 


3. We encounter a block in the dual algorithm applied to Subproblem (1). This 
occurs when a column of Subproblem (1) contains a negative element in row 


k, and a zero in row K;, as illustrated in Figure 2.3. 


Let kg be the column in Subproblem (1) causing the dual block. We now 
define a further contraction, Subproblem (2), and switch to the primal algorithm. 
Subproblem (2) consists of all the rows of Subproblem (1), but only those columns 
of Subproblem (1) with zeros in its target row (kz). Note that the current target 
column in Subproblem (1) (k3) must have a negative element in its “extremal” row 
(k2). As a notational convenience we reverse the sign of this element and record this 
action by labeling the column as —k3. Column —k3 becomes the right-hand side 
column of Subproblem (2), as shown in Figure 2.4. 

The number of columns in Subproblem (2) is reduced by at least one from 
the number of columns in Subproblem (1). We now apply the primal algorithm to 


Subproblem (2) until (which may be after zero or more primal pivots): 





Figure 2.5: Primal Block in Subproblem (2) 


l. vj, ~~, > 0, implying that v4, 4, € 0. À pivot in column з no longer threatens 
to violate primal constraint k4 . We have thus removed the original primal 
block, and we proceed by returning to the original primal problem, and com- 


puting new ratios using columns Kk, and Хз. 


2. v 4, 2 0. A pivot in row k, no longer threatens to violate dual constraint 


ks, and thus we have removed the dual block in Subproblem (1). We revert to 


Subproblem (1) and continue with the dual algorithm. 


J. We encounter a block in the primal algorithm applied to Subproblem (2). We 


illustrate this situation in Figure 2.5. 


Let ks be the row causing the primal block in Subproblem (2). Since the the 
target row in Subproblem (2) (k2) is identically zero for all but the right-hand side 
column (—E3), primal unboundedness cannot occur. Also, the entire ky row in the 
original problem is invariant to pivots in this subproblem. We proceed just as we did 
when faced with this situation in the original primal problem. The current target 
row (k4) in the current problem (Subproblem (2)) becomes the extremal function 
row in a newly defined contraction (Subproblem (3)). All columns of the current 
problem are retained, but only those rows with zeros in the right-hand side column 


of the (—43) current problem are carried forward to the contraction. This again 
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Figure 2.6: Subproblem 3 


assures a strict reduction in the number of rows. The new subproblem is shown in 
Figure 2.6. We then proceed with the dual algorithm, exactly as before. 

The construction of the contracting subproblems through the nested sets of 
zeros in the columns and rows guarantees a monotonic decrease in the sizes of the 
higher-order subproblems. This ensures the ultimate resolution of degeneracy and 
gives us a complete, symmetric and finite algorithm for the solution of LP problems. 

The blocking resolution scheme given here is a constructive algorithm to iden- 
tify strictly improved solutions. The restricted subproblems ultimately yield a pivot 
sequence satisfying all higher-order criteria. Geometrically, we systematically search 
the degenerate subspace for an improved representation. This is in sharp contrast 
to ad hoc “anti-degeneracy” and “anti-cycling” schemes which invoke arbitrary sec- 
ondary mechanisms not at all related to the geometry or mathematics of the problem 
at hand, and consequently admit nuisance degenerate pivots with no constructive 


motivation. 


I. RELATIONSHIP BETWEEN PRIMAL-DUAL ALGORITHM AND 
SIMPLEX METHOD 

We now depart from the presentation of Graves [1965] to discuss the relation- 

ship between the Mutual Primal-Dual Method and the classical Simplex Method. 

After a brief discussion of the similarities, we will explain our reasons for adapting 


the Primal-Dual view. 
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The Simplex Method assumes primal feasibility, and we must identify a start- 
ing basic feasible solution. Because this is seldom practical, the usual approach is to 
formulate and solve a new linear program closely related to the original which has 
the same optimal solution (assuming one exists) and possesses an easily identified 
basic feasible solution. This related problem is derived from (PLP) by augmenting 
it with slack variables, resulting in the following: 

(APLP) min : шу + 08 


s.t.: Ay + 1з=т }т 
Гу >0 }n 
Is>0 }m 


In the classical Simplex view, a basis is a collection of m linearly independent 
columns. Let Ag be such a simplex basis which corresponds to a primal feasi- 
ble solution for (APLP). Suppose we partition the coefficient matrix in a manner 
determined by the basic slack variables, yielding (after possible row and column 


interchanges): 


Ап Ар |1 0 
Ар = 
ад аз де 0). 


where 


_ | An 0 
А = | E | 


and the matrix of nonbasic columns is 


|h Яғ 
Ахь = | "x 


Сә 
СӘ 


The basic variables are yp = (yi, 52) while the nonbasic variables are yyg = 
(31, y2), and the principal part of the Simplex tableau is [Az Ans]. 
Borrowing from the perspective of the Primal-Dual algorithm, we may generate 


the Simplex tableau by partitions: 


Au 0 | 


"ua 22 = n 0 
ІҢ А, h -АлА Lh], 


and the tableau becomes 


ali м Аі 0 I, Aj 
Am] -| tlle 
_ п Ар Ав 
-Ån Á Аз – АЛА Аг |, 


which is precisely the principal part of the Primal-Dual tableau. This is no surprise 
when we note that the basis for the Dual Algorithm, T, is exactly Apg, and the Primal 
and Dual Algorithms share the same tableau. We may interpret the Primal and Dual 
Algorithms as simply two different perspectives of the same tableau, wherein the 
Primal Algorithm a pivot is viewed as exchanging primal constraints and in the 
Dual Algorithm a pivot is viewed as exchanging dual constraints. The classical 
Simplex Method may then be interpreted as solving (PLP) using the 
Dual Algorithm perspective. That the classical Simplex Method is naturally 
interpreted as a dual algorithm comes as a surprise to the conventionally trained. 
However, the consequent mathematical iei is compelling, especially in light of 
the notational simplification and apparent underlying role of Aj)’ . 

There are several reasons for preferring the Primal-Dual Algorithm to the 
Simplex Method. From a computational standpoint, because slack variables are 


carried logically rather than introduced explicitly, we are able to clearly identify the 
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essential information needed to execute the algorithm. The matrix Aj,’ plays a key 
role in the calculation of the tableau, and the entire tableau can be constructed from 
Aj; and original problem data. Since Ajj’ is a submatrix of Ag’, it is smaller and 
requires fewer arithmetic operations to update than does Aj’. 

А second advantage of the Primal-Dual Algorithm lies in the flexibihty it 
offers for specialization to particular problem classes or structures. Indeed, it is the 
special structure and simplicity of the nonnegativity constraints that motivate the 
development of the algorithm in the first place. It is frequently the case that other 
special structures can be identified in classes of (PLP). Examples of such structures 
include simple upper bounds, generalized upper bounds, variable upper bounds, 
pure and generalized network substructures, etc. Such structure may be static in 
that its nature and dimension remains fixed throughout the solution process, or 
the structure may be dynamic in which case its precise nature and/or dimension 
may vary as the problem is solved. Some special structures may be more strongly 
characterized by their column structure and others by their row structure. The 
Primal-Dual Algorithm allows one to effectively exploit virtually any such problem 
structure in a natural manner and greatly simplifies the implementation of such a 
specialization. 

When the linear programming problem appears as a subproblem in a more 
sophisticated solution setting (for example, in a mixed integer programming problem 
or a nonlinear programming problem), the row/column symmetry of the Primal- 
Dual Algorithm is of critical importance in specializing the solution approach. The 
inherent symmetry of the algorithm permits easy adaptation to branch-and-bound 
and cutting-plane approaches to mixed integer programming, to column generation 


settings. as wel] as to primal and dual decomposition techniques. 
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We believe the reason for this flexibility offered by the algorithm lies in its 
more complete mathematical foundation. There is a natural consistency that 
arises from the choice of a vector space having the same dimension as the 
problem variables that is lacking in other approaches. A natural geometric 
interpretation of the solution trajectory follows directly from this development. In- 
cidental issues such as finding an initial basic feasible solution and dealing 
with degeneracy are resolved constructively in this mathematical frame- 
work. Other approaches resort to unnecessarily complicated tangential efforts. 

All the research results reported here can be developed, with some effort, in 
the framework of the classical Simplex Method. However, we choose to present these 
results in the manner of their development - the mutual Primal-Dual view presented 


by Graves. 
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III. IMPLEMENTATION DESIGN OVERVIEW 


A. INTRODUCTION 

We seek to demonstrate the efficiency of row factorization (in particular, using 
dynamic-dimension bases). Accordingly, we will implement the ideas developed here 
and extensively test them within a commercial-quality optimization system: the X- 
System [1975] employs the Graves mutual primal-dual algorithm in a variety of 
large scale optimization applications, including linear, nonlinear, mixed integer and 
decomposed models. The benchmark test suite is drawn from a wide variety of actual 
applications, and our goal is to improve the efficiency of an already well-known and 


highly regarded system. 


B. DESIGN CONSIDERATIONS 

We want to test our ideas by repeatedly solving many medium- to large-size lin- 
ear programming problems (1.e., approximately 8,000-10,000 constraints and 15,000- 
20,000 structural variables). Larger problems are of interest, but for purposes of this 
research, we are limited to a relatively modest computer budget on an IBM 3033/AP 
computer under the VM/CMS operating system using VS Fortran 1.4.1. We wish 
to support the computational enhancements common in commercial mathematical 
programming systems (e.g., bounded variables, ranged constraints, parametric pro- 
gramming, etc.). We require a primal-dual implementation that offers complete 
flexibility in determining solution strategv. In addition, each experimental imple- 
mentation must support all the routine housekeeping of the optimization system 
(e.g., re-inversion, crash starts, relaxations, restrictions, extrinsic and intrinsic enu- 


meration. etc.). Finally. we seek the capability to identify desired row factorization 
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structure within the LP model instance, either by communication from the modeler 


or automatically. 


C. DESIGN TEMPLATE 

To establish a conceptual framework for the evaluation of our algorithms, it 
is useful to outline the important aspects of our implementation and identify the 
crucial steps which most strongly influence performance. Recall from Chapter 2 our 


basic tableau in partitioned form: 


(3) (33) (722) 


(2) m Am A12 Am 
(22) АЛ Án Аі А»2 s Аз Ату А12 JF Á2 ДҮ, (3.1) 
(579 -w Aj} ш – ш А А2 и Агт} 


where we have made the substitution y? — B^!f in the right-hand side column: 


HEN T 
UE at ү ГЫ 11 11 4412 1 
у = В / -» (У1,9) Ec | 0 =. Ю | 


- (wy) = (Ann, 0)’. 
The data requirements of the algorithm are: 


1. Access to the original problem data; 


2. A representation of that part of the current tableau we have chosen to represent 


implicitly, and: 


3. A representation of that part of the current tableau we have chosen to represent 


explicitly. 


We now consider each of these requirements in greater detail. 
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Primal simplex implementations typically are column-oriented and thus re- 
quire column-wise access to the problem data. In the primal-dual method, restric- 
tion to either column-wise or row-wise access alone exacts a serious computational 
penalty. Thus, our design allows both column-wise and row-wise access. 

The design of an efficient large-scale tableau element generation representation 
is remarkably complex. Many subtle software engineering and hardware environment 
issues can have a profound influence on the intricacy and performance of promising 
designs. The design excursions reported here are inexorably influenced by architec- 
ture of the host computer, operating system and implementation language. However, 
the reported design philosophy has been tempered with experience on many other 
computers of widely varied designs. The proposed innovations adapt quite well to 
floating-point pipelines, large cache memories and parallel architectures. 

Because of its fundamental role in the construction of the tableau, we maintain 
an explicit representation of Aj! (we thus refer to Aj as the “explicit kernel”). 
Although any of the various techniques such as LU, LDL! or QR decomposition 
or product form inverse (e.g., Golub and Van Loan [1985] or Magnanti [1976]), are 
suitable for representing A, or Aj, , a difficulty arises in this algorithmic setting. 
While in the primal Simplex method all updates to the basis take the form of rank 
one column exchanges, our setting admits more general updates, which include single 
row exchanges, single row and column exchanges, single row and column deletions, 
and multiple row exchanges of Aj, (multiple column exchanges of A,,). While 
this does not preclude the use of any particular representation, it adds a level of 
complexity not usually encountered in more traditional implementations. 

We also maintain an explicit representation of the right-hand side column and 
the bottom (cost) row. Because of the symmetric nature of the mutual primal-dual 


method, a sensible approach is to allocate a single storage array for both quantities, 
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which taken together are cane’ the “problem rim”. If we generate an explicit rep- 
resentation of the pivot row and pivot column of the tableau at each iteration, then 
we may update the problem rim array using the simple pivot transformation. By 
adopting the convention of labeling the nonbasic constraints as rows 1 through m 
and labeling the basic constraints as rows (m+ 1) through (m+n) (assuming our 
problem instance has m structural constraints and n nonnegativity constraints), the 


problem rim array is partitioned as follows: 


1. The first portion of the array holds values corresponding to the nonbasic con- 
straints in region (2) of the tableau (3.1). Since these are nonnegativity con- 
straints and they are nonbinding (nonbasic), the values in region (2) of the rim 


array are those of the currently (possibly) nonzero variables. 


2. The second portion of the array holds values corresponding to the nonbasic 
structural constraints in region (72), and thus the values in the rim array are 


the current slack or violation in these constraints. 


3. The third region of the rim array, beginning at position m + 1, corresponds 
to basic (binding) structural constraints, and thus the values are those of the 


corresponding (possibly) nonzero dual variables. 


4. The fourth region of the rim array corresponds to basic nonnegativity con- 
straints, and thus the values are those of the corresponding (possibly) nonzero 


dual variables. conventionally called “reduced costs”. 


The rest of the tableau we represent implicitly, by simply recording in a stor- 


age array the current ordering of nonbasic and basic constraints. When additional 


information from thc tableau (for example, a row or column) is required, we con- 


struct it from our representation of Aj,', the current row ordering and the original 


problem data. 
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An overview of the solution process is as follows: 


Identify an initial basic solution. As stated previously, the origin is always 


such a solution, and thus we may always begin with: 


В° = —1 


D? = А, 
or any other suitable basis. 


Check for optimality. If there are currently no primal violations ( a negative 
value in the first or second region of the rim array) and there are no dual 
violations (a positive value in the third or fourth region of the rim array), 


then the current solution is optimal. Otherwise, we proceed to step 3. 
) 


Select either a primal or a dual violation, perform a pivot which makes progress 


towards reducing that violation and return to step 2. 


once we desire to mmaimmtam current information m the rim array by means of 
the pivot transformation updates, we require a representation of the pivot row 
and pivot column at each iteration. Since we explicitly maintain only Ajj! and 
the problem rim, we see that a key computational step in our implementation 


is the generation of tableau rows and columns. 


Recall from Chapter 2 the principal part of the primal tableau (a symmetric 


development from the dual perspective is also possible, and is of course equivalent): 
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(3) (33) 


DP = (1) ; Ап Ay Аз (3.2) 
(ii) \—А А, Ag2 — An Ayy Ai? 
where P = —B™' is the conjugate row basis, 
' Ay, A 
в= 0) | Ап Жа 3.3 
(77) | g o = | -— 


i ғ. - | (3.4) 


D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 
Now consider the generation of column s from (3.2). Rewriting (3.2) in a 


manner that highlights our intentions: 
(2) (39) 


(7) ТЕ АТ A? 


D = (3.5) 


(ii) (Aal Ai) -Aal An A) + Az 


By properly sequencing our computations we will exploit the fact that region 
(31) of a given column is simply a linear combination of region (7) of the same column. 
Assume we want to place the current representation of column s into a work 
array z, which we partition as z7 — (21. 25)7 to correspond to (3.5). Expressed in 
terms of an explicit transformation kernel Aj)’, we first compute region (1) of column 


чав: 


Ый? if s is in (7). 


dm 


t) 


OI 


24, = A Ang): if s is in (07) 


Having done this, we then compute 22 as: 


zg = – Аз 21 if s is in (j), 
or 
22 = – А121 + (А2): if s is in (77). 
Assuming an LU representation of A,, we first compute region (т) of column 
S as: 
DO ate if column s is in (7), 
Or 
LU z, = (Aj2)° {соитип d n D 


Having done this, we then compute 2» as: 


29 = – А121 if column s is in (7), 


OT 


22 = — A223 t (Aga if column s is in (jj). 
Then the current representation of column s is available in 27 = (гу,22)1. 
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E. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN ROW GENERATION 
The computation of row t of the tableau proceeds in a similar manner. We 


now view the principal part of the tableau as: 


(2) (22) 
~m ~n 
DP = (1) Ат (An )А,; 


(i) (AnA (-АлА)А2 + А2 
If we want to place the current representation of row 1 іп а могК аттау 2 


partitioned as 2 = (23,24), we may first compute region (J) of row t as: 


eA if row 1 is in (t), 


OT 


33 = (—An): Al? if row t is in (ii). 


We then compute: 


2 2 (CAO if row t is in (ċ), 


ОГ 


21 = 23(412) + (Аз), if row Ї is in (21). 


Alternately using an LU representation of Aj. 
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AáLU-2e if row t is in (t), 


or 
23LU = (—An)t if row t is in (ii). 
We then compute 24 as: 
2 — 23(Aj2) if row t 1S in is 
Or 
24 = 2з(А\2) + (А»2)‹ if row t is in (ti), 


and the current representation of row t is available in z = (23, 24). 

We see that in each case calculations proceed by first using a representation of 
Aj; to compute a portion of the row or column and then using this initial computa- 
tion and original problem data to compute the remaining part. We will discover that 
our specializations extend this approach by introducing additional tableau partitions 
which allow this computational strategy to be applied on a larger scale. 

As previously mentioned, an important implementation challenge in this algo- 
rithmic setting is the dynamic behavior of A,,. We see from the primal row basis 
(3.3) and nonbasic rows (3.4) that the dimension of A; corresponds to the number 
of basic structural constraints, or, equivalently, to the number of nonbasic nonneg- 


ativity constraints (recall that if a nonnegativity constraint is nonbasic and thus 
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nonbinding, the corresponding variable may possibly be nonzero). Recalling that 
our primal view of a pivot 1s as an exchange of constraints between B and D, we 


see that one of four cases may occur during a pivot: 


1. A structural constraint enters the basis B and a structural constraint leaves 
the basis and enters D. Since the number of basic structural constraints (and 
the number of nonbasic nonnegativity constraints) remains unchanged, the 
dimension of Aj; is unchanged. A pivot of this type involves a row in region 
(j) of B (3.3) and a row in region (22) of D (3.4), and thus it corresponds to 


a pivot coordinate in the location ((72), (7)) of the tableau (3.5). 


2. A nonnegativity constraint enters the basis and a nonnegativity constraint 
leaves the basis. Again, the dimension of A,, remains unchanged. Since this 
pivot involves a row in region (jj) of (3.3) and a row in region (2) of (3.4), the 


corresponding tableau (3.4) pivot coordinate lies in ((2), (j7)). 


3. A structural constraint enters the basis and a nonnegativity constraint leaves 
the basis, and thus the number of basic structural constraints (equivalently, 
the number of nonbasic nonnegativity constraints) increases by one. The di- 
mension of Aj; is increased by one. This corresponds to a pivot coordinate in 


region ((27),(j7)) of the tableau (3.5). 


4. A nonnegativity constraint enters the basis and a structural constraint leaves 
the basis, and thus the dimension of A,, 1s decreased by one. The correspond- 


ing pivot coordinate in (3.5) is ((z), (7)). 


We see that we mav exert some influence on the behavior of the dimension of 
Ay, by our strategy for selecting target violations for primal and dual constraints 


(i.e.. our pricing strategv) and through our tie-breaking rules for choosing pivot 
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row/columns, and that this dynamism is an inherent feature of our algorithm. We 
have already seen the fundamental importance of the kernel (A,;) in our compu- 
tations. Thus, a successful implementation must manage this dynamic behavior 


efficiently and reliably. 


To illustrate in a familiar setting the challenge this offers, consider a LU rep- 


resentation of the kernel Aj): 


A pivot in tableau coordinate ((z), (77)) results in a column exchange in A11. Writing 


(3.6) in the more convenient form: 


alee = U. (Samm 


When column a* replaces the p column of A), to form Aj), we have 


ТІ 
D Ц 
ean = | ’ 
| ^ 
where a — L-!a* has replaced the p'^ column of U. We now must restore the 


upper triangular form of U. We would prefer a method which demonstrates strong 


numerical stability and which also preserves the sparsity of U. Several methods 


сл 
to 


have been proposed (Bartels and Golub [1969], Forrest and Tomlin [1972], Saunders 
[1976] and Reid [1982]). Perhaps the most widely used method is that of Saunders. 
Assume for the moment that Aj, has a block (or *bump") triangular form (we 


discuss how this is done shortly). Then A, appears as shown in Figure 3.1. 





Figure 3.1: Bump Triangular Form of A, 


Each bump consists of (possibly) several columns that extend above the main 
diagonal. These columns are called “spikes”, and the tallest spike within a bump is 
placed in the right-most column, so that a bump appears as shown in Figure 3.2. 

The block triangular form of Aj, results in a LU decomposition which has the 
form shown in Figure 3.3. 

Saunders exploits the form of these LU factors by maintaining a permutation 
of the columns and rows of U so that all the spikes appear on the right-hand side 


without changing their relative ordering. yielding: 





Figure 3.2: Spikes within a Bump 


ай ———————— ---- ---- - e ! 
є ; + 


When the p column of A,, is replaced by a*, the method proceeds as follows: 
р р y 


1. Delete the pt column of U, ù’, and move all the following columns one 


position to the left. 


2. Place a = L'a’ in the right-hand column of Ù. 








Figure 3.3: LU Decomposition of A), 


3. Move the p'^ row of U, up, to the bottom of Ü, and move all rows in between 
up one position. Note that this preserves the triangular form of all but the 
last column of U, and further note that row à, has nonzeros only in columns 


corresponding to R. 


4. Using Gaussian elimination with row interchanges, transforni ù, to a sin- 
gleton column, thereby restoring the upper triangular form of U. Thus, 
U=E,...E,U, where U is the final updated form of U and £,...., £, are the 
elementary transformation matrices corresponding to the Gaussian elimination 


steps (see, e.g.. Murtagh [1981]). 


5. Apply these same Gaussian elimination steps to L^!, forming 


peer EAM 


Typically, L^! is held in product form and thus these transformations are 
simply added to the current list of transformation vectors (“eta” vectors). 

This method has the virtue that new nonzeros can be created only in the 
submatrix F of U, and thus R may be carried in peripheral storage. The numerical 
properties are reported to be quite good. 

The computational burden of continuously maintaining Aj, in block triangular 
form is excessive, and the usual approach is to refresh the structure as part of the 
basis re-inversion routine. Between the re-inversions, the structure is left untended. 
An effective heuristic due to Hellerman and Rarick [1972] is commonly used for this 
purpose. 

A pivot in tableau coordinate ((22), (7)) results in a row exchange in Aj. 


Writing (3.6) as 


im EE (3.8) 


When row a, replaces the q'^ row of Aj, to form Aj, we have 


where row Â = a,U^! replaces the q'^ row of L. The desired structure of L may be 
restored by methods symmetric to those developed for the column exchange case, 


again forming : 
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"m == LU . 
A pivot in tableau coordinate ((22),(j7)) causes the dimension of Aj; to in- 
crease by one through the addition of a row and a column. It is convenient to add 


the new row to the bottom of Aj, and the new column to its right. If Aj, is the 


new kernel, then 


А К 
Ay = | A” | 


a, а; к 


(3.9) 


where a,,a* and а, are original matrix coefficients. 


The desired updated decomposition is of the form: 


Ai ak = L 0 U u* 
ab as E ND UH. 


which requires that 


and 


k 
ат Е = T x Tz ne 


Setting uz. = 1. we have: 


К 
pe — ак ЖР l.u . 


The final pivot coordinate, ((2), (7)) results in the dimension of Aj; decreasing 


by one. If Aj, is the current transformation kernel and Ay, results after the pivot, 


without loss of generality we have: 


an а.к а; E dum 0 Urr М, 
Ay, = | ас танде | = | FL | li U, | | (3.10) 


where (ak, a,) is the leaving row, (а, к, a*)? the leaving column and a,, the 
pivot element. Using an analysis similar to Rosen [1960], we note that if the square 


matrix C is nonsingular and is partitioned as: 


С © 
di: A , 


where C is square and nonsingular, C, is square and Cy = Cy — C3Cy'C) is non- 


singular, then 


Cc Qr ОТОС ето а 
- (СОС б 


Applying this result to (3.10), we discover that 


Ашу ш а 


ОГ 





^ 1 
Ay = L,U; + | 
а; к 


ES (3.11) 


П the term (==) a^a. should be the zero matrix, our new decomposition is at 


а, к 


hand. Unfortunately. (=) aa, need not be zero, but we may guarantee 1t to 


be so by performing a preliminary column exchange update to Аң. We replace 


5 


СУ) 


Шет: 
column (ar, a*) with column (1, 0)? and comipute the updated factors exactly as 
we did in the column exchange case considered earlier. This results in a modified 


transformation kernel Ан with factors 


А өн 1 n Е m 0 err ü, 
ы 0 А} И l; Lı 0 U; 


Then the second term of (3.11) is 


1 
(=) 0-4, =0 


and thus the updated LU factors are at hand: 


a == eC, . 


We now see the advantages and disadvantages of the Graves mutual primal- 
dual method. It has the advantage that, although the transformation kernel can 
be as large as m, the number of structural constraints in the problem instance, its 
dimension 1s actually equal to the current number of binding structural constraints 
(or equivalently, to the number of potentially nonzero primal variables). At our 
initial basic solution, the kernel dimension is zero. If the maximum kernel dimen- 
sion during the course of the solution trajectory is much smaller than m, we enjoy 
significant computational advantages in storage, update and computation. For a 
great many LP model instances, this is precisely the case. 

The dynamic nature of the transformation kernel clearly complicates update 
procedures. While the usual primal Simplex method requires only the column ex- 


change update. Graves’ method requires four update cases. 
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In our implementation, we will seek methods for handling the transformation 
kernel that enjoy the advantages while mitigating the disadvantages. 

It is this general strategy for row and column generation we develop and en- 
hance in our specializations. By identifying a special structure in the problem data, 
we will introduce additional linear dependencies in the rows and columns of the 
tableau which will further simplify their generation and will also further reduce the 


dimension of the transformation kernel. 
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IV. FACTORIZATION 


A. INTRODUCTION 

Each of the algorithms we present in this research can be viewed as a special- 
ization of a general approach to large-scale linear programming developed by Graves 
and McBride [1976], which they call the “factorization approach”. It is based on 
distinguishing special rows and columns in a way that allows large parts of the ba- 
sic tableau to be represented implicitly, to be generated easily from the remaining 
explicit parts only when actually required by the algorithm. 

Although sometimes used interchangeably in the literature, we recognize a 
distinction between partitioning methods and factorization methods. In the former, 
a formal distinction is made between substructures which appear in the model in- 
stance, usually constraints or variables. An approach for solving the problem is then 
developed which exploits the substructures. For example, Dantzig-Wolfe decompo- 
sition is such an approach which partitions constraints, while in Benders decompo- 
sition the variables are partitioned. Such an approach may be applied statically, in 
which case the desired partitions are identified at the start of the solution process 
and remain fixed throughout, or it may be dynamic, in which case the partitioning 
may be adjusted as the solution proceeds. 

In contrast, we consider factorization methods to be those in which a formal 
distinction is made between substructures in the LP basis (and thus in the basic 
tableau). The algorithm is then specialized to exploit this special substructure. 
Thus, we view the GUB algorithm of Dantzig and Van Slyke [1967] and the mul- 


ticommodity network flow algorithm of Hartman and Lasdon [1970] as examples 


01 


of factorization algorithms. Factorization methods may also be static or dynamic, 
depending on whether the dimension of the factored substructures in the LP basis 
may vary during the solution process. 

Our research is concerned with dynamic row factorization algorithms. We. 
develop the general dynamic factorization approach here and discuss the important 
aspects of the algorithm. In subsequent chapters, we specialize this general approach 
to each of several special LP row structures. This general development closely 


parallels that of Graves and McBride [1976]. 


B. THE FACTORED TABLEAU 


The problem to be considered 1s: 


(FLP) min : wy 


Sd шс T } explicit constraints 
Ғи<% }factored constraints 
—Iy <0 )nonnegativity constraints , 


where y is a n vector of decision variables, w is a n vector of cost coefficients, E 
is an m by n matrix of constraint coefficients for "explicit" constraints with right- 
hand side m vector r, F is a p by n matrix of constraint coefficients for “factored” 
constraints with right-hand side p vector b, and —7 is the negative of the n by n 
identity matrix. Іп this general development, we refer to the F-type constraints 
as “factored” only to distinguish them from the “explicit” E-type constraints, and 
assume nothing about their structure. Not until our specializations in later chapters 
will we impose special structure on F, and the structures we will consider permit 
the representation of the F-type constraints without the inversion of a matrix. We 
will see that this approach is centered around handling the part of the basis cor- 


responding to the E-type constraints explicitly while factoring the portion of the 
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basis corresponding to the F-type constraints. The notation is chosen to suggest 
this idea. 

As we saw in Chapter 2, in the mutual primal-dual method the primal al- 
gorithm and the dual algorithm share the same tableau, and thus the tableau for 
(FLP) may be derived from either perspective. We present only the derivation from 
the primal viewpoint here. 

Recall that a basis for the primal algorithm consists of n linearly independent 
rows from the constraint coefficient matrix when it is assumed to include both 
structural (explicit and factored) and nonnegativity constraints. Assume that the 
current row basis consists of | rows from E, k rows from F and (n — k — l) rows 
from —J. Using the notation of Chapter 2 temporarily, we have: 


l+k п-і-К 
= к= 


pe [An An | M+R 
0 — І \п- 1-Е, 
where [A A12] includes all basic structural rows, from both E and F. 
We will ultimately be interested in isolating the effect of each type of struc- 
tural constraint algebraically in the factored tableau, and thus we require greater 
resolution in our factored basis. Introducing obvious notation, we have: 


К l п-і-К 
е aie жағы 


En Еһ Ез |}! 


[Au Ai) — 
Е; Fiz Fis } k › 
where 
Ea En 
EIE 
Ан E Fa 
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From the development in Chapter 2, we recognize A4; as the explicit kernel, 
and thus Aj,’ exists. It then follows that it is always possible to identify among the 
columns of [Fi; Fi2] a nonsingular submatrix Fy, of dimension k, since otherwise 
the rank of [Fi Fiz] is at most (k — 1) and thus the rank of B is at most (n — 1). 
We will later see that one of the important implementation challenges is the task of 
efficiently managing the structure and nonsingularity of Fi. 


The full factored row basis 1s then: 


k l n—-l-k 
~m~ N m 


(7) E E Ey }1 
B= G3) | Pa Fa № | м. 
ОТТО 0 -І )n-i-k 


Introducing the notation: 


^ A^ E 
Ay = Ез 29 ЕЕЕ} 


ZN = 
Ay = E3 — En Fi! Fis 


2 


its inverse 1s: 


ЕРА (1+ Fy FA En) Fa’ Fy) (Fis — Fi Aj A2) 
B! = Ат – Ат Еу Рп. Ay Ai 
0 0 —[ 
Grouping the coefficients of the nonbasic constraints and applying the same 


column ordering vields: 
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(i) [-1 0 0 Е 
ШЕ (а) Ға ЕЁ, Рз |р —k 
(4)| 0 -I 0 )! 


(9) (Ед En Ез }т—1 
We will explain the ordering of the rows of D shortly. 


(4.2) 


The principal part of the factored tableau is DP, where P = —B™ is the 


conjugate row basis. Introducing the additional notation: 


Fy = Fo. — Fx Fy’ Fi2 
ee ee ee 
лил ЛЫ Ee 
А; = з ЛУД. 


the principal part of the factored tableau 15: 


(7) (27) (722) 
ЖЫ еле М 
(Q) -Fi Från (I+Fp Frå Ena) fn Fo (F -— FAnn A2) 
mp (22) ПІ д! (А Ет ES uA E — TAI 
(iii) Ату -А En Fy! Ар Ау 
(iv) \ TARA аи — En) Fū’ boy = PAR 
(4.3) 
Partitioning w — (wj.w3,ws,w4), 7? — (rj,r;)! and b = (b,b) and 


introducing the notation: 


А = 
U aE w — wiii Fiz 
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-1 
Ша = и Ет Ез 


103 = 

bac b,- Fe 
ry = rı — En Fū’ b 
T) = r2 — En Fb , 


the complete factored tableau 15: 


A (1+ Fū Fn An En) Fa Ет (Fia E Fi Ату Аз) Fi (b im Ату?) 


-Fn AÑ (Fn An En — Fa) Fi Fog — Fr Avy Án b, — Рз Аут 

Ат а | Ann Án fi (4.4) 
-АлАр (Ал Ар Еп — Ex )Fy’ An = Ал Ат А12 no An An Tı 
-шАр (із АТ Еу = w,) FR! os Ф; Атү А12 ui Fb гр Фу Ату Тү 


We see from (4.4) that with knowledge of the current factorization, we can 
construct the entire tableau from Fj;!, Aj! and the original problem data. The 
dimension of Fj,’ is equal to the number of F-type constraints that are currently 
basic, and thus can be at most p. The dimension of Ajj! is equal to the number of 
E-type constraints that are currently basic. Hence, its dimension cannot exceed m. 
We call Ajj! the explicit transformation kernel and Fū’ the factored transformation 
kernel. 

We have seen in Chapter 3 that the origin is always a convenient initial basic 
solution for the mutual primal-dual method, and the same is true for the factoriza- 


tion approach. An initial basic solution is always: 


GG 


Starting from this solution, we may view the algorithm as progressing by exchanging 
rows of F and E from D with nonnegativity constraints from B. We will find it useful 
to associate with each F-type constraint in B a unique nonnegativity constraint in 
D. Borrowing the terminology of Dantzig and Van Slyke [1967], we will call each such 
unique nonnegativity constraint (and its corresponding variable) the "key" variable 
for the associated F-type constraint. The algebraic structure arising from this basic 
F-type constraint/nonbasic *key" variable association allows us to represent parts 
of the tableau implicitly rather that explicitly. To distinguish these *key" variables 
in D from others, we place them in region (2) of D. This 1s the reason for the row 


ordering mentioned earlier. 


C. BENEFITS OF FACTORIZATION 
Graves and McBride [1976] report three principal benefits of the factorization 


approach: 
1. Good starting bases; 
2. Reduction in memorv requirements; 
3. Reduction in work per pivot. 


Although the origin is always available as an initial basic solution, Graves and 


McBride [1976] suggest that a better one can frequently be found with only small 
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computational effort using this approach. They form the Lagrangian Dual of (FLP) 


with respect to the E-type constraints, we obtain: 


(LFLP) min: wy + ХЕу-т) 
5.0.7 ЖШ 5 


—Iy <0 

where \ > 0 is a m vector of dual variables. If À equals A*, the optimal dual variables 
for (FLP), then an optimal solution of (FLP) is among the optimal solutions of 
(LFLP). They speculate that if Å is a “good” estimate of \*, then using À in (LFLP) 
and solving for y should yield a “good” starting point ý for (FLP). If no such estimate 
À is available, then À 2 0 may be used. Depending upon the structure of the F-type 
constraints, the resulting problem may either be much easier to solve than (FLP) and 
thus may be solved optimally, or heuristics may be used to find a computationally 
inexpensive “good” solution. All factorization algorithms can be started by solving 
(LFLP) first. 

Standard simplex techniques require a basis of dimension (m + p) to solve 
(FLP). The factorization approach requires the explicit transformation kernel Aga 
whose dimension is at most т, апа the factored kernel Fi;, whose dimension is 
at most p. Thus. it is clear that we expect a considerable reduction in memory 
requirements using the factorization approach. When we specialize the approach 
for models in which the F-type constraints have special structure (e.g., GUB, pure 
network or generalized network), we will see that memory requirements can be 
further reduced. 

While it 1s difficult to measure the computational effort per pivot when product 
form inverse or basis decomposition techniques are used, there are indications that 


the factorization approach can vield substantial benefits. Both analytic (McBride 
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[1972]) and empirical (Tomlin [1972]) studies indicate that when the F-type con- 
straints are GUB constraints (we will cover this in detail in Chapter 6), the factoriza- 
tion approach results in per-pivot computational benefits when the GUB constraints 
comprise at least 3096 of the total number of structural constraints. With some sim- 
plifying assumptions, a similar analysis can be developed for the case where the 
F-type constraints are pure network rows. 

We will consider each of three different factored row structures in this work. 
The complexity of the factored row structure determines the computational difficulty 
of representing and updating the factored kernel Fj;. The simplest row structure 
we consider, in which the factored rows are generalized upper bounds, allows F}; to 
be interpreted as a simple permutation matrix. The second factorization, in which 
the F-type constraints are pure network rows, 1s a more complex structure which 
subsumes generalized upper bounds as a special case. In this case, Ff}; may be 
interpreted as representing a partial ordering defined over a subset of the rows of F. 
The fina] factorization, in which the rows of F are generalized network rows, is the 
most general we consider in this work and subsumes both the other factorizations. 
Неге, Г тау Бе interpreted as a linear transformation in which a near partial 


ordering exists among a subset of the rows of F. 
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V. GENERAL IMPLEMENTATION TOOLS 


A. INTRODUCTION 

We now provide an overview of our implementation of the factorization method. 
We describe the representation and update of the three important algorithm compo- 
nents: the basic tableau, the factored kernel and the explicit transformation kernel. 
The fundamental update steps are described in terms of functions which operate 
on each of the algorithm components. Where the details of the update actions are 
common to all three factorizations treated in this research, the details are presented 


in this chapter; otherwise, they appear ir subsequent chapters. 


B. THE FACTORED TABLEAU 
As before, we are interested in the problem: 


min: wy 
y 


s.t.: | Ey € r) explicit constraints 
Fy < 6} factored constraints 
—Iy< 0} nonnegativity constraints 


Recall the algebraic representation of an arbitrary primal row basis: 


(7) En En Es 
(77) Ға ЕЕ; Р» (5.1) 


ыл ат 21 


ШІ -- 


with the corresponding nonbasic rows: 


ч) "Mm 
p- (i) | Fa Fa Fa 
(iiij o -I o0 
G AEn En DA 


The conjugate row basis inverse is: 


(5.3) 


Fu Fun Еп РРА ЕР -Fa P Fa FoAn (Eis - En Fn P) 
P=-B'= -An Aii En Fn’ -А (Ез — ЕР Аз) 
0 0 1 


and the principal part of the factored tableau is: 


(i) амы Е ЫЕ А ЫБ БС EcL ERE ENATUCE E uESEE) 
(i) | -Frå 4 FuFü Fus РА ЕР КАРД! Fo- Fa Fa Po + Fin Fa’ FAn (2 - Eu Fi’ Fis) 
ZPJ EA -ҒрАТЦҒЕа- ҒаҒа Ға) 
(iit) Ат -Ân En Fi’ Ån (En - En Ез) 
ne MEE LU EaR od EnA EuFa' - Enf Ға ЕЕЕ АКЕ ВЕРИ) 
| ЕЕ РА ЕТЕ See he EnA En = Eika Pi) 


The entire principal part of the factored tableau can be constructed with 
knowledge of the row and column ordering, Ёп! (ог Ё.) and Áil. Our discussion 
of the implementation centers around the representation and update of these three 
components. 

We restrict our attention to the principal part of the tableau because our 
implementation maintains a current representation of the right-hand side column 
and bottom (cost) row (collectively referred to as the "problem rim") at all times. 
These quantities are always immediately available and need not be computed upon 
demand. This is in contrast to (5.4), where only Az} and Fy; are always available 
and any other quantity of interest (for example, a complete row or column) must 


be computed. 


C. ALGORITHM OVERVIEW 


In broad terms, the algorithm proceeds as follows: 


1. Establish an initial basic solution. 


2. Stop if the current solution is optimal. Optimality exists when the right-hand 


side column is nonnegative and the bottom (cost) row is nonpositive. 


3. Based on a pricing scheme and a ratio test, select and generate a row and 


column of (5.4). 
4. Stop if this row and column reveal infeasibility or unboundedness. 


5. Transform the problem rim. the tableau representation, F,; and Aj}, and 


réetummeto step 2 


After developing the data structures and necessary update operations, we will 


expand upon this overview. 


--1 
М 


D. IMPLEMENTATION CONVENTIONS 

We will first discuss the implementation of each of the three algorithm com- 
ponents individually, and then provide a detailed presentation of the complete al- 
gorithm. Our implementation is coded in the FORTRAN language, and thus the 
following discussion is colored by the character of that language. 


1. The Tableau 


The mutual primal-dual method fosters a symmetric view of the linear 
programming problem and the solution process. This symmetry is reflected in the 
indexing conventions required in our implementation. Assuming now that (5.1) 
consists of m structural constraints (both factored and explicit) and n structural 
variables, we require the structural constraints to be indexed from 1 to m, and the 
variables from (m+ 1) to (m+n). While our statement of the problem specifies 
the nonnegativity constraints explicitly, we will of course deal with them implicitly. 
Note, however, that we have two equivalent interpretations of the indices (m + 
1) through (m+n). We can view them as the indices of the structural problem 
variables or as the indices of the nonnegativity constraints corresponding to each 
of the structural variables. We will use both interpretations in our development. 
Observe also that we require no explicit representation of logical (slack, surplus or 
artificial) variables. 

The indexing of the original problem data exists independently of the 
solution process and is thus referred to as the “extrinsic” indexing. The factoring of 
the constraints to form B and D represents a second indexing of the problem which 
defines the current solution, independent of the extrinsic indexing. We refer to this 
indexing system as “intrinsic”, and it partially defines the current representation of 
B and D. Our convention for intrinsic indexing has the rows of D labeled from 1 


to m and the rows of B labeled from (m+ 1) to (m+n). 
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The original problem data is stored in super-sparse form, (in super-sparse 
form, each real value is stored once, and each nonzero coefficient is represented by a 
row index, a column index, and a pointer to a real value) with unique nonzero real 
values stored once in an array of type DOUBLE PRECISION, and referenced by the 
FORTRAN equivalent of a pointer from each constraint matrix coefficient. These 
pointers to the real values are singly-linked by both row (with associated column 
indices) and by column (with associated row indices). This allows symmetric row- 
wise and column-wise access to the data to the performed efficiently, which is an 
important virtue in a primal-dual design. 

To specify the structure of (5.4) we require a representation of the cur- 
rent mapping between the extrinsic indexing of the problem data and the intrinsic 
indexing of the current tableau. We represent this mapping by two arrays of type 
INTEGER, each dimensioned from 1 to (m+n), which are used as pointers. We 
denote these arrays MINT() and MEXT(). MINT() represents the mapping from 
intrinsic to extrinsic indices, and MEXT() represents the inverse mapping. Thus, 
if IINT is an intrinsic index of a row or column of (5.4) and IEXT is the extrinsic 
index of the constraint or variable currently mapped to position IINT of the tableau, 


then: 


MINT(IINT) = JEXT 
MEXT(IEXT) И 
MINT(MEX TEXT) TETEXT 
MEAT MATTIA TD “ІЛЕ 
To complete the representation of (5.4) we require knowledge of the fac- 
torization which is represented algebraically by regions (7), (ii). (iii). etc. We 
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handle this with the use of several variables of type INTEGER which are used as 
pointers to record the ending intrinsic index of the various tableau regions. Table 


(5.1) lists the regions and the associated pointer names. 


TABLE 5.1: Indices of Tableau Regions 


REGION BEGINNING INDEX ENDING INDEX 


(1) 1 MKC 
(11) MKC + 1 MFR 
(111) МЕК +1 MXR 
(iv) MXR + 1 M 
j) eee NXR 
(jj) NXR +1 NFR 
(33) МЕК +1 ШЕ 


The arrays MINT() and MEXT() and the pointers MKC, MFR, MXR, 
NXR and NFR comprise the data structures necessary to represent the principal 
part of the tableau. To complete the representation of the tableau we require a rep- 
resentation of the right-hand side column, the bottom row and the current value of 
the extremal function. The right-hand side column and bottom row are represented 
by an array of type DOUBLE PRECISION, dimensioned from 1 to (m+n), referred 
to as RIM(). Positions 1 through m correspond to the right-hand side column and 
positions (m+ 1) through (m +m) correspond to the bottom row. Note that these 
are intrinsic indices. The current value of the extremal function is held in a scalar 
variable of type DOUBLE PRECISION called OPT. 

We now discuss the operations necessary to maintain and update the 
tableau representation. To recognize an optimal solution we simply scan the RIM() 
аттау. When the values in positions 1 through m are nonnegative and those in 


positions (m+ 1) through (m+n) are nonpositive, the current solution is optimal. 


=] 
Qt 


Violations of either sign discipline may serve as potential primal algorithm target 
rows or dual algorithm target columns, respectively. 

Having selected a target row or column, we must generate the correspond- 
ing row or column of the tableau, compute ratios with respect to the right-hand side 
column or bottom row, and thus identify the index of a corresponding column or 
row. We then generate that column or row of the tableau. Because the details of row 
and column generation differ according to the factorization, we defer their discussion 
until the appropriate chapter. To represent the update operations, however, we de- 
fine the abstract functions Generate Row(IPRI) and Generate. Column(IPCI). 
These functions accept as arguments the intrinsic row (IPRI) or column (IPCI) 
index of the current tableau and return the current representation of that row or 
column. We require a working array of type DOUBLE PRECISION to hold these 
current representations. We define ROWCOL/() to be such an array, dimensioned 
from 1 to (m 4- n), where positions 1 through m correspond to column IPCI of (5.4) 
and positions (m 4- 1) through (m -- n) correspond to row IPRI. ROWCOL() is 
indexed intrinsically in conformance with our convention for (5.4). 

Recall that our interpretation of a pivot from the primal algorithm view- 
point is as an exchange of rows between B and D. The row indices of D correspond 
to the row indices of (5.4), while the row indices of B correspond to column indices 
of (5.4). Thus, a pivot requiring the exchange of a row of B and a row of D requires 
the exchange of a row index of (5.4) with a column index of (5.4). Using MINT() 
and MEXT(), such an exchange requires just a few assignment statements. For 
example. if the pivot row IPRI corresponds to extrinsic index EPRI and the pivot 


column IPCI corresponds to extrinsic index EPCI, then prior to a pivot we have: 


MINT(IPRI) - EPRI 
MEXT(EPRI) - IPRI 
MINT(IPCI) = ЕРСІ 


МЕХТ(ЕРСІ) = ІРСІ 


and after the pivotal exchange we have: 


МІМТ(ІРЕІ) = EPCI 
MEXT(EPCI) = IPRI 
MINT(IPCI) = EPRI 


MEXT(EPRI) = IPCI. 


We call this type of index exchange a “primary exchange”. For abstrac- 
tion purposes we define a function Primary_Index_Exchange (IPRI,IPCI) which, 
given the current tableau representation and the intrinsic pivot row index IPRI and 
the intrinsic pivot column index IPCI, performs the appropriate update of MINT() 
and MEXTY(). Every pivot requires a primary exchange. 

In addition to the partitioning of constraints into B and D, we are re- 
quired to maintain the factorizations ((1), (12), etc.) within each. Having performed 
the primary exchange. it is possible that additional exchanges within B and/or D 
are required to restore the proper factorization. For example, if the primary ex- 
change removes an F-type constraint from D (region (?:)) and places it in region (7) 
of B, an additional exchange is necessary to move the row from region (J) to region 
(23). We refer to such actions as ^secondary exchanges". Secondary exchanges in D 
correspond to row exchanges in (5.4) and secondary exchanges in B to column ex- 


changes in (5.4). Thus, we define the functions Row Index Exchange(IRI1.IRI2) 


те) 
=! 


and Column.Index.Exchange(ICI1,ICI2) which exchange intrinsic rows IRI1 and 
IRI2 or intrinsic columns ICI1 and ICI2, respectively. Since the indexing of the 
RIM() array corresponds to the tableau indexing, we assume that any secondary 
exchanges performed on the tableau are also performed on RIM(). 

Finally, our factorization requires that the factored kernel, F31, remain 
nonsingular. It is possible that the primary and secondary exchanges will destroy 
the nonsingularity of Fj,;. When this occurs, additional row exchanges between 
regions (i) and (iit) of D will be necessary to restore the required nonsingularity of 
Fi. We call such exchanges "tertiary exchanges". Again, we assume that tertiary 
exchanges performed on the tableau are also performed on RIM(). 

The factorizations of B and D are dynamic in nature, meaning that 
a particular pivot may cause the dimension of one or more regions of B and/or 
D to increase or decrease by one row. Such dimension changes are handled by 
simply adjusting the values of the factorization pointers. We define the functions 
Increment(XXX) and Decrement(XXX), where *XXX" is MKC, MFR, MXR, 
NXR or NFR, to increment or decrement, respectively, the appropriate pointer 


value. 
In summary, the tableau data structures are: 


MINT() 
MEXT() 
RIM() 
ROWCOL() 


MKC, MFR, MXR, NFR, NXR, 


~] 
On 


and the necessary update operations are: 


Generate Row(IPRTI) 

Generate Column(/IPCT) 
Primary Index Exchange(I PRI,IPCI) 
Row Index Exchange(1RIl,IRI2) 
Column.Index.Exchange(ICI1,ICI2) 
Increment(X X X), and 


Decrement(X X X ). 


2. The Explicit Transformation Kernel 

As can be seen from (5.4), the rows of Aj] correspond to region (iii) of 
the tableau (nonbasic nonkey variables) and the columns to region (j) (basic explicit 
constraints). The dimension of Ajj) may vary as the algorithm progresses, and the 
mechanism for these changes is the addition and deletion of rows and columns. 
We thus require data structures that permit the efficient insertion, deletion and 
exchange of rows and columns of Aj. Further, the number of nonzero elements 
in Aj} varies from pivot to pivot, independently of dimension changes. We store 
these nonzero elements in a stack, referring to them via pointers which are stored 
as an orthogonally linked list, doubly linked by row and doubly linked by column. 
This allows the efficient update of Aj! and accommodates the dynamics, at some 
expense in storage. Contrary to the tableau representation, we maintain the explicit 
transformation kernel representation indexed extrinsically rather than intrinsically. 


This allows permutations within region (j) and/or within region (277) of the tableau 


without requiring corresponding permutations in the columns and/or rows of the 
Aj; representation. 

We again define a suite of functions which operate on our representation of 
Az} and describe its update. Add_Row(IRE) and Add_Column(IJE) append ex- 
trinsic row IRE and extrinsic column IJE, respectively, to our representation of Aj}. 
We will ensure that the current representation of each is in work array ROWCOL() 
prior to executing these functions. Columns will always be appended to the right- 
hand side of Aj} and rows will be appended to the bottom. Delete_Row(IRE) and 
Delete_Column(IJE) delete extrinsic row IRE and IJE, respectively, from the rep- 
resentation of Aj). Replace-Row(IRE1,IRE2) causes the overlay of extrinsic row 
IRE1 by extrinsic row IRE2. The representation of IRE2 will have been previously 
placed in ROWCOL/(). Finally, Update_Explicit_Transformation_Kernel per- 
forms a pivot update of the representation of Ajj, using the current representation 
of the pivot row and pivot column in ROWCOL(). 

3. The Factored Kernel 

The data structures and update actions for the factored kernel vary ac- 
cording to the factorization and their discussion will be deferred. Here we define the 
necessary abstract update functions. 

Is Factored Kernel.Singular tests the current representation of Fy, 
for singularity and returns the appropriate Boolean value. It will be used in certain 
pivot cases to determine which of several courses of action should be taken. 

All tertiary exchanges take the iur of row exchanges between regions 
(7) and (277) of (5.4). Note that region (7) consists of the nonnegativity constraints 
of key variables. By our alternate interpretation of the nonnegativity constraint 
indices we may consider these to be key variable indices. Since the columns of F4 


are key variables. we interpret region (7) as containing the indices of the columns 
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of F\,. An exchange of rows between (7) and (272) of (5.4) is then interpreted as an 
exchange of columns of Fj. 

Before such a column exchange for Fi, (and thus a row exchange in (5.4)) 
can be made, we must frequently identify the index of the column to be removed from 
Fy, (which is an index in region (:) of (5.4)) and the index of the column which will 
replace it (which is an index in region (777) of (5.4)). We thus define Find_Column 
-to.Remove(IOUT), which selects from among the indices of region (2) of (5.4) 
the intrinsic index IOUT of the column to be removed from Ё. Similarly, Find 
-Column_to_Add(IIN) selects from among the indices of region (tti) of (5.4) the 
intrinsic index IIN of the column to be be added to F};. Finally, Update_Factored 


-Kernel updates the data structures used to represent Fi. 


E. COMPLETE ALGORITHM DESCRIPTION 


We are now in a position to fully describe our implementation of the 
algorithm. We do so by expanding the discussion of each step given previously in 


section C. 


1. Initialize the algorithm using the origin as the initial basic solution. Then 


LU 
^- (E 
ve = (F) 


2. Stop if the current solution is optimal. The current solution is optimal if 


(RIM(IR) » 0for1 € IR € M and RIM(JC) <OforM+1<JCSM+N. 


3. (a) Select a violated primal or dual constraint as the target row or column, 


respectively. For purposes of exposition. assume the current solution is 


el 


C 





4. 


5. 


primal feasible and we are executing the primal algorithm. Then the 
target row is the bottom row, and we select the intrinsic index IPCI of a 
column which will allow us to make a gain in the value of the extremal 


function (i.e., RIM(IPCI) » 0). 


Generate.Column(IPCI). 


ENT 
4 


Calculate ratios by computing RIM(IR)/ ROWCOL(IR) for those 1 € 


— 
eo 
— 


IR< Mwith ROWCOL(IR) > 0. If ROWCOL(UIR) < 0 forall JR,1 < 
ІВ < M stop, since the problem is primal unbounded. Otherwise, select 
IPRI to be the intrinsic index yielding the smallest such ratio. Break 
ties in accordance with the following hierarchy: region (222) first, then 
region (ii), then (i) and finally region (iv). Within a region, break ties 


by selecting IPRI to be the first such index encountered. 


(d) Generate Row(IPRI). 


If, contrary to our assumption in 2.a, the current solution 1s not primal feasible, 
the target row would be some row IPRI of (5.4), 1 € IP RI € M, rather than 
the bottom row. We would proceed by next executing Generate Row(IPRI). 


If ROWCOL(JC) > 0 for all Af +1 < JC < M+N, we would conclude that 


the problem is primal infeasible. 
(a) Primary_Index_Exchange(IPRI,IPCI). 
(b) Pivot update RIM() and OPT using ROWCOL(). 


(c) Perform the necessary secondary and tertiary exchanges, as shown in Ta- 
ble 2. Note that some tertiary update actions are conditional, depending 
upon the singularity of Fj;. We use a notation similar to the FORTRAN 


BLOCK IF statement to indicate the conditional actions. 
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(d) Update.Factored. Kernel. 
(e) Update. Explicit Transformation Kernel. 


(f) Go to Step 2. 


We now give a more detailed explanation of the secondary and tertiary 


exchanges listed in Table 5.2, discussing each pivot coordinate individually. 


1. Pivotal coordinate ((1), (j)). EPCI, the extrinsic index of an E-type constraint 
located in position IPCI of B, is exchanged with EPRI, the extrinsic index 
of a nonnegativity constraint located in position IPRI of D. The initial pivot 
exchange places EPCI in position IPRI of (2) and EPRI in position IPCI of 
(7), and thus secondary row and column exchanges are necessary to restore the 
structure of the tableau. Beginning with the column structure, we exchange 
column EPRI in position IPCI of region (7) with the extrinsic index located in 
position NXR of region (j77), and then decrement the value of NXR. This has 
the effect of shifting the starting column index of region (jj) one position to 
the left. Now extrinsic index EPRI is in position (NXR+1) of the tableau and 
thus is still misplaced. We then exchange EPRI in position (NXR+1) with the 
extrinsic index located in position NFR. This places EPRI contiguously with 
region (JJJ), so by decrementing NFR we restore the column factorization. 
The effect has been to decrease the dimension of (J) by one, increase the 
dimension of (jjj) by one and to shift the entire (JJ) region one position 
to the left. In the row structure, EPCI in position IPRI of region (i) must 
eventually be moved to region (iv). However, prior to this pivot EPRI was a 
key variable, located in region (7). Its removal has disturbed the structure and 
thus the nonsingularity of F\,;. We must therefore identify a column currently 


in region (?i?) which can be used to restore the nonsingularity of Fjj.. We 
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thus invoke Find_Column_to_Add(IIN), which identifies the intrinsic index 
IIN of such a column. Suppose the extrinsic index of the column in position 
IIN is KX. We exchange EPCI in position IPRI of region (7) with KX in 


position IIN of region (zii). We now have extrinsic index EPCI in position 


IIN of (iii) and we must move it to region (tv). We thus exchange EPCI | 


in position IIN in region (iii) with the extrinsic index in position NXR of 
(iii), and then decrement NXR. This completes the restructuring of the row 
factorization with the effect of decreasing the dimension of region (272) while 
increasing that of region (iv). The effect of these exchanges has been to reduce 
the dimension of Aj} through the deletion of column EPCI and row KX. To 
maintain the proper structure of Aj, we therefore Delete Row(KX) and 


Delete. Column(EPCI). 


. Pivotal coordinate ((7), (77)). EPCI, the extrinsic index of a basic F-type con- 
straint located in position IPCI of B, is exchanged with EPRI, the extrinsic in- 
dex of a nonbasic nonnegativity constraint located in position IPRI of D. Since 
both EPCI and EPRI are misplaced after the primary exchange, secondary row 
and column exchanges are needed. In the columns, we exchange EPRI in po- 
sition IPCI of region (77) with the extrinsic index located in position NFR of 
region (7j). Decrementing NFR then completes the column exchanges, leaving 
the dimension of region (jj) reduced by one and that of (jjj) increased by 
one. In the rows. we exchange EPCI in position IPRI of region (7) with the ex- 
trinsic index in position MKC of region (i). Decrementing MKC then restores 
the row factored structure of the tableau. It remains to be seen, however, 
whether what remains of F}; after the removal of row EPCI and column EPRI 


is still nonsingular. We test F1; by invoking Is.Factored. Kernel. Singular. 


54 


If the response is FALSE, then no tertiary exchanges are required. If, how- 
ever, the response is TRUE, we must identify a column exchange that we will 
restore the nonsingularity of Fj;. To identify the two columns which will be in- 
volved in the exchange, we invoke Find. Column.to.Remove(IOUT) which 
returns the intrinsic index IOUT (with associated extrinsic index NX, say) of 
a column to be removed from F1;. We then call Find_Column_to_Add(IIN) 
which identifies the intrinsic index IIN (with associated extrinsic index KX, 
say) of a column in region (iii) which can be used to replace IOUT. We now in- 
voke Row Index. Exchange(IIN,IOUT), which completes the update of the 
row factorization. This row exchange results in the replacement of row KX of 
by row NX of the tableau. We currently have no representation of row KX of 
the tableau, so we compute it by invoking Generate Row(IOUT). We may 


then restore the structure of Åj} by invoking Replace. Row(KX,NX). 


. Pivotal coordinate ((1),(737)). EPCI, the extrinsic index of a basic nonneg- 
ativity constraint located in position IPCI of B, is exchanged with extrinsic 
index EPRI located in position IPRI of D. The primary exchange properly 
places EPRI in region (777), and thus no secondary column exchanges are 
needed. Likewise, EPCI is properly placed in position IPRI of region (1), so 
no secondary row exchanges are required. Prior to the pivot, EPRI was desig- 
nated a key column and its removal disturbed the structure of Fj;. Since EPCI 
is itself a column, it is possible that it may be a replacement for EPRI. To 
determine if this is the case, we invoke Is_Factored_Kernel_Singular. If the 
response is FALSE, the direct exchange of column EPCI for column EPRI in 
is justified and the tableau update is complete. If the response is TRUE, how- 


ever. a tertiary row exchange is required. EPCI, currently located in position 


IPRI of region (:), will be the column we remove from F1, and its replacement 
is found by invoking Find.Column.to. Add(IIN) which returns the intrinsic 
index IIN (associated with extrinsic index KX, say) located in region (227). We 
perform a tertiary exchange by invoking Row_Index_Exchange(IIN,IPRI). 
This completes the row factorization update. This tertiary exchange results 
in the replacement of row KX of Áil by row EPCI of the tableau. We must 
therefore perform the corresponding row exchange in Áil also. To perform 
such an exchange, we require the current representation of row IPRI of the 
tableau. This is precisely the pivot row and is already available in the work 


array ROWCOL(). Thus, we may invoke Replace. Row(KX,EPCI) directly. 


. Pivotal coordinate ((27),(7)). EPCI, the extrinsic index of a basic E-type 
constraint located in position IPCI of B, is exchanged with EPRI, the ex- 
trinsic index of a nonbasic F-type constraint located in position IPRI of D. 
After the primary exchange, both EPCI (in (4:2)) and EPRI (in (7)) are mis- 
placed, and thus both secondary row and column exchanges are necessary. 
We exchange EPRI in position IPCI of region (7) with the extrinsic index 
in position NXR of region (7) and then decrement NXR. This has the ef- 
fect of reducing the dimension of region (7) while increasing that of region 
(77). Notice that we have added row EPRI to the structure of ХР}, апа 
thus an additional column must at some point be identified. Several row 
exchanges are needed. First, we exchange row EPCI in position IPRI of 
region (77) with the extrinsic index located in position (MKC+1) of region 
(ii) and then decrement MKC. This restores the structure of region (27) by 
reducing its dimension by one row. We now must execute a series of ex- 


changes that ultimately places EPCI in region (?v) and also adds a column 
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to region (z) to restore the nonsingularity of Fj;. To do this, we first iden- 
tify the column to be added to F), by invoking Find_Column_to_Add(IIN), 
which identifies the intrinsic index IIN in region (iii) (associated with an 
extrinsic index KX, say) which may be added to F}. Notice that the dimen- 
sion of Fi; is increased by one through the addition of row EPRI and col- 
umn KX. We then perform Row_Index_Exchange(IIN,MKC), which places 
KX in position MKC of region (i) (correctly) and places EPCI in position 
IIN of region (272) (incorrectly). To complete the tertiary row exchanges, we 
invoke Row_Index_Exchange(IIN,MXR), followed by Decrement(MXR). 
This properly places EPCI in region (iv) by increasing the dimension of re- 
gion (iv) by one row while decreasing that of region (iii). Finally, notice that 
the dimension of Aj} has been reduced by one through the removal of column 
EPCI and row KX. To update the representation of Aj} accordingly, we invoke 


Delete Row(KX) followed by Delete Column(EPCI). 


. Pivotal coordinate ((72),(77)). EPCI, the extrinsic index of a basic F-type 
constraint located in position IPCI of B, is exchanged with EPRI, the ex- 
trinsic index of a nonbasic F-type constraint located in position IPRI of D. 
The primary exchange places both EPRI and EPCI in their proper positions, 
so no secondary row or column exchanges are needed. However, the struc- 
ture of has been altered through the addition of row EPRI and the deletion 
of row EPC]. To determine if the structure and nonsingularity of Fy; has re- 
mained intact, we invoke Is. Factored Kernel.Singular. If the response is 
FALSE, the tableau update is complete. If the response is TRUE, we must 
identify a column exchange that will restore nonsingulanty. Thus, we invoke 


Find_Column_to_Remove(JOUT) which returns the intrinsic index ОСТ 


Ср) 
-a 


(associated with the extrinsic index NX, say) of a column located in region 
(i) which may be removed from Fj,. We invoke Find. Column.to.- Add(IIN) 
which returns the intrinsic index (associated with an extrinsic index KX, say) 
of the column IIN located in region (tii) which may be added to Fy,. We then 
perform a tertiary exchange by invoking Row.Index Exchange(IIN,IOUT). 
This last action results in the replacement of row IIN of Aj} by row IOUT of 
the tableau. We have no current representation of row IOUT, and to compute 
one we invoke Generate. Row(IOUT). We then update the representation of 


Aj} by invoking Replace_-Row(KX,NX). 


. Pivotal coordinate ((22), (jj7)). EPCI, the extrinsic index of a basic nonnega- 
tivity constraint located in position IPCI of B, is exchanged with EPRI, the 
extrinsic index of a nonbasic F’-type constraint located in position IPRI of D. 
The primary exchange misplaces both EPRI and EPCI and thus both sec- 
ondary row and column exchanges are necessary. We exchange EPRI, located 
in position IPCI of region (jjj) with the extrinsic index located in position 
NFR. By then incrementing NFR, we restore the column factorization struc- 
ture by increasing the dimension of region (77) while reducing that of region 
(777). Notice that the addition of row EPRI to the structure of Fy; indicates 
that an a corresponding column must also be located. To restore the structure 
of region (11), we exchange EP CIin position IPRI of region (7) with the extrin- 
sic index in position MXC of region (2). Incrementing MKC restores the struc- 
ture of region (11) by reducing its dimension by one row. Since EPCI is a col- 
umn which we have now placed in position MKC, it is possible that F; is now 
nonsingular. To determine this, we invoke Is_Factored_Kernel_Singular. 


If the response is FALSE, the tableau is complete. If, on the other hand, 


58 


the response is TRUE, we must perform a tertiary column exchange to re- 
store the nonsingularity of Fjj. We will be removing EPCI, currently located 
in position MKC, from fF); and replacing it with the column identified by 
Find.Column.to.Add(IIN), which is extrinsic index KX located in posi- 
tion IIN of region (iii). We then perform Row.Index Exchange(IIN,MKC), 
which completes the tertiary exchange. This exchange implies a row exchange 
іп Агу also. We are required to replace row KX of Aj}! with the current 
representation of row EPCI of the tableau. Since IPRI is the pivot row, the 
current representation of EPCI is readily available in ROWCOL(). We may 


then complete the row exchange by invoking Replace Row(KX,EPCI). 


. Pivotal coordinate ((222),(7)). EPCI, the extrinsic index of a basic E-type 
constraint located in position IPCI of B, is exchanged with EPRI, the extrin- 
sic index of a nonbasic nonnegativity constraint located in position IPRI of 
D. The primary exchange misplaces both EPRI and EPCI, and thus both 
secondary row and column exchanges are necessary. We exchange EPRI, in 
position IPCI of region (7), with the extrinsic index in position NXR of re- 
gion (j). We then exchange EPRI, now in position NXR of region (7), with 
the extrinsic index located in position NFR of region (77). We complete the 
column update by decrementing both NXR and NFR, which has the effect of 
reducing the dimension of region (J) by one column, shifting region (JJ) one 
column to the left in the tableau and increasing the dimension of region (777) 
by one column. For the row update, we first exchange EPCI, in position IPRI 
of region (izi), with the extrinsic index located in position MXR of region 
(722). By then decrementing MXR we restore the row factorization by increas- 


ing the dimension of region (iv) by one row while decreasing that of region 
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(iii). The effect on Ajj has been to reduce its dimension by one through the 
deletion of column EPCI and row EPRI. The corresponding update to the 
representation of Ajj! requires that we invoke Delete Row(EPRI) followed 


by Delete_Column(EPCI). 


. Pivotal coordinate ((722),(j7)). EPCI, the extrinsic index of a basic F-type 
constraint located in position IPCI of B, is exchanged with EPRI, the extrin- 
sic index of a nonbasic nonnegativity constraint located in position IPRI of 
D. The primary exchange misplaces both EPRI and EPCI, and thus both 
secondary row and column exchanges are necessary. We exchange EPRI, cur- 
rently located in position IPCI of region (77), with the extrinsic index currently 
located in position NFR of region (7j). By then decrementing NFR, we re- 
store the column factorization by decreasing the dimension of region (j7) while 
increasing that of region (777). Note that the removal of EPCI from the struc- 
ture of Fi, implies that the dimension of Fi, decreases by one, and thus a 
corresponding column must be removed from Fj. To locate this column, we 
invoke Find.Column.to. Remove(IOUT), which returns the intrinsic index 
IOUT (associated with the extrinsic index NX, say) of a column located in 
region (1) which may be removed from Fj; while allowing the remaining struc- 
ture to be nonsingular. We exchange EPCI, currently located in position IPRI 
of region (277), with NX, currently located in position IOUT of region (1). This 
restores the row structure of region (122), but EPCI is misplaced in region (2). 
Therefore, we exchange EPC] in position IOUT of region (7) with the extrin- 
sic index located in position MKC of region (7). By then decrementing MKC, 
we restore the structure of region (7) by decreasing its dimension and that 


of region (27) by increasing its dimension. The exchange of NX in IIN with 


90 


10. 


E. 


EPCI in IPRI implies a row exchange in Aj. To perform this exchange in 
the representation of Aj, we must replace row EPRI of Az} with the current 
representation of row NX of the tableau. Since this current representation is 
not available, we invoke Generate_Row(IIN) to compute it. We then update 


our representation of Aj} by invoking Replace_Row(EPRI,NX). 


. Pivotal coordinate ((222),(jj7)). EPCI, the extrinsic index of a basic nonneg- 


ativity constraint currently located in position IPCI of B, is exchanged with 
EPRI, the extrinsic index of a nonbasic nonnegativity constraint currently lo- 
cated in position IPRI of D. The primary exchange properly places EPRI and 
EPCI, so no secondary exchanges are necessary. Since F} is unaffected, no 


additional action is required. 


Pivotal coordinate ((iv),(J))- EPCI, the extrinsic index of a basic E-type 
constraint currently located in position IPCI of B, is exchanged with EPRI, the 
extrinsic index of a nonbasic E-type constraint currently located in position 
IPRI of D. The primary exchange properly places EPRI and EPCI, so no 
secondary exchanges are necessary. Pj; is unaffected, so no additional action 


is required. 


Pivotal coordinate ((iv),(j7)). EPCI, the extrinsic index of a basic F-type 
constraint currently located in position IPCI of B, is exchanged with EPRI, 
the extrinsic index of a nonbasic E-type constraint currently located in po- 
sition IPRI of D. The primary exchange misplaces both EPRI and EPCI, 
so both secondary row and column exchanges are necessary. EPRI, cur- 
rently in position IPCI of region (7j). is exchanged with the extrinsic in- 
dex located at position (NXR+1) of region (77). The column structure is 


restored by incrementing NXR. increasing the dimension of region (7) while 


О 


decreasing that of region (jj). Note that since the dimension of region (7) 
increases, the dimension of will increase also. Also, the removal of EPCI from 
Fi; implies that the dimension of Fj; will decrease. Thus, a column of Fi 
must be identified for removal as well. To find such a column, we invoke 
Find.Column.to.Remove(IOUT), which returns the intrinsic index IOUT 
(associated with the extrinsic index NX, say) of a column whose removal from 
F; allows the remaining structure of F}; to be nonsingular. We perform a ter- 
tiary exchange by exchanging NX, currently located in position IOUT of region 
(1), with the extrinsic index currently located in position MKC of region (+). 
The structure of region (?) 1s then restored by decrementing MKC, reducing 
the dimension of region (1) while increasing that of region (1:1). To properly po- 
sition EPCI and restore the structure of region (12), we exchange NX, currently 
located in position (MKC41) of region (22), with EPCI, currently located in 
position IPRI of region (7v). This leaves only NX misplaced. We thus exchange 
NX, currently located in position IPRI of region (iv), with the extrinsic index 
located in position MXR of region (27). The row update is completed by incre- 
menting MXR, increasing the dimension of region (222) while decreasing that 
of region (iv). Since the dimension of Aj} has been increased, we must up- 
date our representation. We are required to add the current representation of 
column EPRI, which is available as the pivot column in ROWCOL(). We also 
require the current representation of row NX of the tableau, which is not cur- 
rently available. We thus invoke Generate Row(IOUT) With these two rep- 
resentations, we may then invoke Add. Row(NX) and Add. Column(EPCI) 


to complete the update of the representation of Aj). 


12. Pivotal coordinate (iv), 20) EPCI, the extrinsic index of a basic nonneg- 
ativity constraint currently located in position IPCI of B, is exchanged with 
EPRI, the extrinsic index of a nonbasic E-type constraint currently located in 
position IPRI of D. The primary exchange misplaces both EPRI and EPCI. 
To restore the column factorization, we first exchange EPRI, located in posi- 
tion IPCI of region (777), with the extrinsic index located in position (NFR+1) 
of region (777). We then exchange EPCI, now in position (NFR+1) of region 
(777), with the extrinsic index located in position (NXR+1). We complete 
the column update by incrementing both NXR and NFR, which has the effect 
of increasing the dimension of region (j), shifting region (JJ) one position to 
the right and decreasing the dimension of region (jjj). Since the dimension of 
region (j) has been increased, the dimension of Aj; will increase as well. To 
restore the row factorization, we exchange EPCI, currently in position IPRI 
of region (iv), with the extrinsic index located in position (MXR+1) of re- 
gion (iv). The row update is completed by incrementing MXR, which has 
the effect of increasing the dimension of region (277) while decreasing that of 
region (?v). The dimension of Aj} has been increased through the addition of 
column EPCI and row EPRI. Since IPCI is the pivot column and IPRI is the 
pivot row, each is already available in ROWCOL(). Thus we simply invoke 
Add.Row(EPRI) followed by Add.Column(EPCI) to complete our update 


of the representation of Aj,’. 
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Secondary and Tertiary Tableau Exchanges 


TABLE 5.2 
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VI. FACTORIZATION OF GENERALIZED 
UPPER BOUND ROWS 


A. INTRODUCTION 


We now develop the first of three specializations of the factorization approach. 


We are still interested in the problem: 


(GUB) min : wy 
st.: Ғу<%8 F pbynm 
Ey<r E mbyn 
-Іу<0-І nbyn 
where F, E,—I,w,b and r are as before. We now require that the F-type constraints 
are generalized upper bound (GUB) constraints. Define 5;,2 — 1,...,p to be pair- 
wise disjoint subsets of the set N — (1,...,n) and further define $, = N — Ü 5; 
ҮШ Бегетаріу). Тһеп 5:Г|5; - 10г0 =: <р, 0 <17 <р, 19 г) 
Ü S; = N. GUB constraints are of the form 


= 


Аы сл 5-7 1,2250) 00 = +1 (6.1) 
265. 
The sets S;, 2 = 0,...,p are called GUB sets and : is the GUB set index. S, 
identifies those variables that have no coefficient in any GUB constraint. The E- 
type constraints are of arbitrary form. 
A specialization of the Simplex Method to handle GUB constraints was first 
developed by Dantzig and Van Slyke [1967]. They introduce a specialization of 
(GUB) with é; = l and strict equality constraints. Their algorithm requires a 


working basis of dimension (m + 1). McBride [1972] develops a specialization of 


the mutual primal-dual method to solve a second variant of (GUB) with 6,, = +1 
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and strict equality constraints. А dynamic working basis whose dimension cannot 
exceed m is required. A computational analysis by McBride [1972] predicts that the 
performance of this algorithm should exceed that of the Simplex Method when the 
proportion of GUB constraints in the model exceeds approximately 29%, and his 
empirical results tend to confirm this analysis. 

We extend this work by developing a specialization of the factorization ap- 
proach which does not require all GUB constraints to be binding. We will see that 
there are computational efficiencies to be gained by this approach. We first present 
the factored tableau for (GUB). We then discuss the important computational is- 


sues. 


B. THE FACTORED TABLEAU 
Recall from Chapter 4 the primal row basis at an arbitrary point in the solution 


process 1s: 


k ! n-l—k 
жз тет a= 


(7) Еп Рл Eis М 
а ОЗГЕ }К (oo) 
А0 0 Іш 
We saw in Chapter 4 that it is always possible to identify among the columns of 
B a set of k columns such that the submatrix Fj, is nonsingular. Since the rows 
of [Fii Fiz Fis] correspond to GUB constraints, their form is similar to (6.1), with 
column permutations accounting for the difference. Each column of {F}; Fiz Fial is 
either zero, or a signed unit vector, and thus by row permutations of these F-type 


constraints, F;; may be placed into the form of a k by k signed identity matrix Aj. 


The resulting (GUB) row basis is: 
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(7) ЕГЕП ЕБ ГІ 
= (77) Ав Fiz Fa } k (6.3) 
Л Оо = pce) oe 


The corresponding nonbasic constraint rows are then: 


k l n-l-k 
Же om 


ie OO 
п Із Жа тен 

ҒА -І 017 

(1) Ena Еп Bo Jie} im —! 


Recall that we wish to observe a one-to-one association between the variables 


(6.4) 


whose nonnegativity constraints (which appear in region (2) of D) are nonbasic 
(nonbinding) and the F-type constraints in region (7j) of B. This association is 
exactly analogous to the idea of “key” variables proposed by Dantzig and Van Slyke 
[1967]. We call the variables whose nonnegativity constraints appear in region (2) 
of D “key” variables, and there is one such variable for each currently basic F- 
type constraint appearing in region (j7) of B. The variables whose nonnegativity 
constraints appear in region (222) of D are then called *non-key". 

Using (6.3) and (6.4) and recalling the expression for the explicit transforma- 


tion kernel: 


ЕЕ БаР) = (Ei = ВИ)" 


9, 


the principal part of the factored tableau is: 


(i) -AuFuádi An t Aun Fidi Enn OnFis — OnFiaAn (£3 — En dufis) 
(11) -Fnán ҒаАт Енде Ға- ҒаЛТҚЕз- Enn Fay) 
(in Ат! -А Елди Ant (Eis — En. Our Fis) 


. a _ 
(ы) | -(Ез- EnduFiw)An —Endu t (En Enn idn Enðn Ens- ҒадҺаРҒіз Е EnA (En < Enôn Fo) 
£EnAu Fan (Ei 7 Endu F3) 


(6.5) 


C. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 


To examine the actions required by Generate Column(LC), suppose we 


want to generate column LC of the principal part of the tableau (6.5) and place 


the results in the vector z7 = (2), 22, 23, 24)1 which is partitioned to conform to 


(6.5). Rewriting (6.5) in a more convenient manner, we have: 


95 


(i) -AuFn [A] -AuFu |-Ai £u] Ап -AuFa [Ai (Es — Еи да Ға) + даР 
(ii) -Fn [Äri] -Fn |-Ā7 Enn] -Fn (АЕ - Елдар) + Рә 
(iii) An -АШЕцда Ay (Eis — En AuFis) 

(iv) | -En [ÂȘ}] - En [-An Falin) -Ев |-А, Еи ди] -En |Ai (£s - £i Fio) 


-En -Anfal ÂR En ôn) + ^u] -Bn [-óu Fut Án Us - БадаРҒз))% AnFis] + En 


(6.6) 


To highlight those computations involving terms consisting entirely of zeros and plus 
and minus ones, we introduce the following notation: a general matrix product is 
denoted by *", as in Aj - Eja. А "simple" matrix product (i.e., one in which one 
of the terms consists entirely of zeros, plus ones and minus ones) is denoted by o, as 
in Ej © Fi;. The calculations may then proceed as follows (note that e”© is a unit 


vector with a plus one in position LC): 
l. Calculate region (277): 


(a) If column LC is in (7), then: 


(c) If column LC is in (377), then: 


23 = P (ғаз) 2 - Е19 24119 (Ғаз) (| 


Notice that (F13)"€ is either null or it is a signed unit vector, say in row 


К in which case Еһ © Ан o Оа = +(Е\,)*. 
2. Calculate region (1): 


(a) If column LC is in (7), then: 


2] — — 1 o Е 9 23 


(b) If column LC is in (77), then: 


Zn AMOR Oo e (Аһ) 


(c) If column LC is in (777), then: 


21) = —Ay,O F429 23 + Ay Oo (Fae 


Notice that the computation — 71; o z4 involves only additions and sub- 


tractions, and that each case in region (:) differs only by an addition 


suffix. 


3. Calculate region (77): 
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(a) If column LC is in (7), then: 


22 = – Е) © 23 


(b) If column LC is in (7j), then: 


Z} = — F5 9 23 


(c) If column LC is in (757), then: 


z2 = — Fn © z3 + (Роз) 
4. Calculate region (iv): 


(a) If column LC is in (7), then: 


СА —En: - Еҙ2: 23 


(b) If column LC is in (7j). then: 


24 = – Е · 21 - Еҙ2: 23 


IBI colum LC is in (777), then: 


z4 = —Ёзү + 2) — EQ za (E53) 
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Then column LC of (6.6) is available as z^ — (z1, za, 23, z4)7. Notice that 
floating point multiplications and/or divisions are necessary only for computation of 
regions (2:1) and (iv); also note that regions (2), (22) and (zv) are linear combinations 
of region (iti) and one another, and that each term in these regions differs only by 


an additive suffix. 


D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN ROW GENERATION 
Now we examine the actions of Generate_Row(LR). Suppose we desire to 
generate row LR and place the results of the computation in the vector 2 — (25, 26, 27) 


which is partitioned to conform to (6.5). Rewriting (6.5) in a convenient form: 


U) (22) (222) 
~ ~ ~ 


(i) ZA o - {[-3n Fn A7] En =- 1] ân [-àu£udi]] Eis» [- (In Fai] £u 9 1) 54] Fs 
(ii) EL - {!-Е»Аү | Еи} 5и [- Fai] Eis 9 [- ([- F2] En] 8i] + А 
(iii) A - ([Az2] £u) 5 [An] Eae [in] £i) 5) 5s 

(ww) (nni - ExlÁ - (Ennii - EM] En + En} Su I En An Fi - БЛАЦЕ, 


+ [- Es an fis -— Ез) | Еһ + Ex} An] Ру ~ Еу 
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Then we may proceed as follows (note that erg is a unit vector with a plus one in 


position LR): 
1. Calculate region (7): 


(a) If row LR is in (1), then: 


£; 2 [-(Au)rn o ix] o A] 


(b) If row LR is in (22), then: 


25 = (—Fo)no9 du 


(c) If row LR is 1n (122), then: 


25 = erro Ay)’ = (Aj) )iR 


(d) If row LR is in (tv), then: 


E (Елін o Ano Fi- (Ев) zen 


Notice that the computation of the term (E»)rg o Aq o Fi; — (Ez )LR 


involves only additions and subtractions. 
2. Calculate region (77): 


(a) If row LR is in (7), then: 


26 59 б; “Еп Кеім)о Аһ 
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(b) If row LR is in (ii), then: 


Že = (—25- Еһ) oan 


(с) If row LR is in (272), then: 


ж-(-25: Еһ) oA, 


(d) If row LR is in (iv), then: 


26 = (—25 Ea — (Ез) oA 
3. Calculate region (j77): 


(a) If row LR is in (2), then: 


27 = 25 · Буз + 260 Рз 


(b) If row LR is in (zz), then: 


ғы 


7 = 25 · Ез + 260 Ёз + (Fo3)LR 


(c) If row LR is in (222), then: 


С cero rae 
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(URR issin- (dv), then: 


2; = 25. Ёз + 26 © Fis + (Ez3)LR. 


The row LR of (6.7) is available as Z — (25, 2, 27). Note that floating point 
multiplications and/or divisions are avoided for terms involving Fi, F2; and F2;, 
that regions (Jj) and (jjj) are linear combinations of region (j), and that many 
cases differ only by additive infix terms. 

This computational scheme extends the approach of the mutual primal-dual 
method by introducing additional linear dependencies into the tableau and exploits 
those dependencies to improve efficiency. Empirical evidence suggests that this 


approach also improves numerical stability (McBride [1989]). 


E. DATA STRUCTURES 


The essential information contained in the factored kernel when the factored 
rows are GUB constraints 1s simply the unique mapping between each basic factored 
constraint and its corresponding "key" nonbasic variable. We require a represen- 
tation of this mapping which is compact, can be efficiently accessed and is easily 
updated. 

One possibility is to maintain the intrinsic ordering of the tableau in such a way 
that the factored row/key variable relationship can be derived implicitly. Region 
(77) of the tableau corresponds to the basic factored constraints of B. Region (:) 
of the tableau corresponds to the nonnegativity constraints associated with the key 
variables. If we maintain the ordering of (7j) and (1) so that the k'^ position of (1) 
contains the index of the kev variable associated with the factored constraint whose 


index is in the kt? position of (jj), we have an implicit representation of Fi. 
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Two difficulties arise with such an approach which lead us to choose an alter- 
nate representation. First, such an ordering complicates the tableau update scheme 
presented in Chapter 5 and potentially increases the computational expense of all 
tableau updates which affect F,,. Second, several update cases require the iden- 
tification of GUB set membership for variables which are not currently assigned 
to region (7) and thus are not currently part of the F\, structure. This may be 
done by accessing the original problem data by column and scanning that column 
for the index of the GUB row, if any, in which the column has a nonzero element. 
Since the original problem data is stored in a super-sparse form, this scheme re- 
quires indirect computer memory addressing. We prefer a method which requires 
less computational expense. 

We thus introduce an additional data structure to manage the structure of 
Fi. KEY(IJ) is an array of type INTEGER, dimensioned from 1 to (m+n), 
which for each basic factored constraint records the extrinsic index of the associated 
key variable, and for each nonbasic key variable records the extrinsic index of the 
corresponding basic factored constraint. Additionally, for each variable which is not 
key (i.e., either non-key variables or basic variables), KEY(IJ) records the extrinsic 
index of the factored constraint for which it may serve as a key variable. If we 
interpret the extrinsic index of the factored constraints as GUB set indices, then 
for every variable KEY(IJ) contains the GUB set index. Note that since GUB 
set membership is fixed for a given variable, all column indices KEY(IJ) can be 
initialized once and need never be updated. Only the factored constraint indices 


require updating. 


F. FACTORED KERNEL UPDATE ACTIONS 


Recall from Chapter 5 the details of the factored kernel update actions 
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Is Factored. Kernel. Singular, 
Find_Column_To_Remove(IF KC), 
Find_Column_To-Add(JF KC), and 


Update_Factored_Kernel, 


are factorization-specific. We now discuss these details for the GUB factorization. 

The row-permuted diagonal structure of Fi; (^11) implies that the question of 
singularity arising from rank-one updates can be resolved strictly on the basis of local 
information. That is, the exchange, deletion or addition of a row leads immediately 
to the identification of a unique GUB set. The search for a corresponding column 
(variable) to complete the exchange, deletion or addition is then limited to those 
sharing this GUB set membership. The simplicity of this structure greatly enhances 
the efficiency of the necessary update actions. 

For example, consider pivot update case ((22),(j7)), in which basic factored 
constraint EPCI, located in tableau position IPCI is removed from the row basis 
and nonbasic factored constraint EPRI, located in tableau position IPRI replaces it. 
Because EPCI was basic, there is a key variable, say EPCIKV, located in position 
IPCIKV of region (1) of the tableau. Because GUB sets are, by definition, pair- 
wise disjoint, EPCIKV may not serve as a key variable for the new basic factored 
constraint EPRI. Thus, a replacement variable among those in region (tti) of the 
tableau must be found. Such a replacement must exist, for otherwise B is singular. 
To locate such a variable, we scan the variables in region (227), searching for one 
which is a member of GUB set EPRI. The test is simply this: for each variable JC 
in region (277), if 

KES JC) bP RE 
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then accept JC as the key variable for EPRI and proceed with the secondary tableau 
update exchange. Otherwise, continue the scan. 

Thus, the action Find.Column.- To Add(EKV) scans the indices of region 
(iii), performing the GUB set membership test on each index. The first such index 
for which the test is true is chosen as EKV. 

The action Find.Column. To Remove(EKV ) is required when a constraint 
ERI is removed from F,,;. EKV is then the key variable corresponding to ERI and 


is determined by: 
EKV - KEY(ERI). 


Is. Factored Kernel.Singular is required when the constraint EPRI is to be 
added to Fi, and the pivot column index EPCI corresponds to a variable. Thus, 
EPCI is a possible key variable for EPRI and is a convenient candidate. It may be 
accepted as the key variable for EPRI if and only if F1; remains nonsingular after the 
addition of both EPRI and EPCI, and this is the case if EPCI is a member of GUB 
set EPRI. Thus, Is_Factored_Kernel_Singular simply compares KEY(EPCI) and 
EPRI. If they are unequal, Fj; is singular and the result “true” is returned. Other- 
wise, the result is “false”. 

Finally, we consider Update_Factored_Kernel. Since the GUB set mem- 
bership of a variable does not change, only constraint index updates of KEY(IR) 
are required. For example, if the old key variable index for constraint EPRI is 


EPRIKVO and the new index is EPRIKVN, the required update is: 


KEY(EPRI) = EPRIKVN. 


VII. FACTORIZATION OF PURE NETWORK 
ROWS 


A. INTRODUCTION 


We remain interested in the problem: 


(PNSC) min [o wy 

s.t.: Fy<b ; F pbyn 

Ey<r ; E mbym 

—[y €0 ; -Inbym 
where F, E, —I, ш, b and r are as before. We now require that each column 
of F consist of zero, one or two nonzero elements. If a column contains a single 
nonzero element within the rows of F, it may be either a plus one or a minus 
one. If a column contains two nonzero elements within the rows of F, one must 
be a plus one and the other a minus one. F-type constraints may also be in- 
terpreted as defining the node-arc incidence matrix of a directed graph. Define 
І = {0,...,p},7 = {1,... n}, and N = {n;,i € Z — (0]) to bea set of nodes and 
A fa, 6.27 | аб-(п,п;,і61,76 1) to be a set of arcs with ordered pairs 
of nodes (tail,head) as elements indexed by k. Note that we interpret node 0 as a 
null node, and thus an arc incident to node 0 is viewed as a single-ended arc. Then 
a graph is defined as G = {N’, A}. The node-arc incidence matrix of € is a matrix 
F consisting of p rows (one for each node in A’) and n columns (one for each arc in 


A) with elements: 


+1 if a* =(n;,n;) for some n; € A 
іл 11а = (n; n;i) for some n; E N 
0 otherwise 


109 


Thus, if arc a* 2 (ni,n;) is represented by column k of F, then 


and each column of F is either all zeros, contains exactly one nonzero element (which 
may be either plus one or minus one) or contains exactly two nonzero elements (a 
plus one and a minus one). 

A specialization of (PNSC) that has been widely studied arises when the F- 
type constraints are taken to be equalities and the E-type constraints are vacuous, 


resulting in: 


(PNE) min 10 wy 
өр: Ке ЛЫ 
—Iy<0 -I nbyn 
Very efficient specializations of the primal Simplex Method have been developed 
and implemented to solve (PNE) (e.g., Srinivasan and Thompson [1973], Glover, et. 
al. [(1974], Bradley, Brown and Graves [1977]). These algorithms exploit some well 
known properties of F (we assume that one redundant row has been removed from 


E 


1. Every primal Simplex basis B of F (consisting of (p — 1) linearly independent 


columns of F) can be triangulated by row and column permutations, 
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2. Every such basis B is itself the node-arc incidence matrix of a subgraph of G, 


and this subgraph is a rooted spanning tree, and 


3. F is totally unimodular, implying that if bj and b; are integer (p — 1)-vectors 
апа =; апа =з are (p — 1)-vectors of unknowns, then for every primal Simplex 


basis B of F the solutions of Bz, — b апа тЇ В = B are also integer. 


Property (1) allows very efficient execution of Simplex iterations, property (2) mo- 
tivates the use of special data structures which allow efficient storage and update of 
simplex bases and property (3) allows all computations to be performed in integer 
arithmetic (assuming w and b are integer). 

Several contributions have been made to solving variations of (PNSC) when 
the E-type constraints are not vacuous. These approaches have their roots in work 
by Kaul [1965] and Bennett [1966]. Chen and Saigal considered a version of (PNSC) 
with all equality constraints. Hartman and Lasdon [1972] considered a specialization 
of (PNSC) to the multicommodity capacitated transshipment problem (MCTP), in 
which the E-type constraints are generalized upper bound (GUB) constraints and 
the F-type constraints decouple into independent pure network subproblems. In 
each of these treatments, the resulting algorithm is a specialization of the primal 
Simplex Method. McBride [1972] considered (PNSC) in the context of the mutual 
primal-dual method, with all constraints assumed to be equalities. Each of these 
approaches allows a basis representation which may be factored into two compo- 
nents; an explicit part of dimension m by m and a factored part of dimension p by 
р. 

Our inequality formulation of (PNSC) allows a dynamic basis representation 
where, just as in the GUB specialization, both the dimension of the factored part 


and the dimension of the explicit part may vary from one iteration to the next. 
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B. THE FACTORED TABLEAU 


Recall from Chapter 3 the primal row basis at an arbitrary point in the solution 


process 1s: 


(j) f En En En 
B= шы ae (7.1) 
ТШ” 


The corresponding nonbasic constraint rows are then: 


() Ше та 
Se) E oS 
ТІ а ш 


(1v) Ea En Ea 
and the principal part of the factored tableau is: 


(72) (0022) 
(5) M —Á 


(1) ETAT p ЕТТЕРІ ЕТЕП ЕП Ёз — Fa! Fadu (En 7 En Fi Ға) 
(й) | -Fazi + FaFa Padu  FandnEunFn —f£afn o Ға- Pu Fo) Fis + Fafa (Ез = Би) 


SENEC Fd EE EPa Es- Erh Fa) 

(iii) An a MES ADM £13 — Eu Fa Fis) 

(iv) | -En ATi + Pafo bala Endi Eum - Ез Е E;5 t En Fi Fin (Eis 7 En Fi Fia) 
zr pupa T =F, Pa Fa Боле Боа 


(7.3) 


The F-type constraints represent the node-arc incidence matrix of a graph. For 
notational convenience, let us define the corresponding directed graph explicitly as 
G = {N, A} where G,N and A are as previously defined. Note that while we choose 


to interpret all columns of (PNSC) as arcs in G, any number of these columns may be 
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null with respect to the F-type constraints and thus may be interpreted graphically 
as null arcs. 

Since the rows of [Fj1. Fi? Fj3] in the primal row basis B form a subset of the 
rows of F,[Fj Fi; Fia] may itself be interpreted as the node-arc incidence matrix 
of a graph GB = {NB, AB}. NB is a subset of NV, but in general Ag ¢ A. To 
see this, define Np = N — Np and consider an arc a? = (n;,n;) where n; € Ng 
and n; € Np. Then a* “spans” the partitioning of NV into Vg and Np. While a? 
is a doubleton arc with respect to G, it is a singleton arc with respect to Gg. We 
call such an arc a “dynamic singleton”, and we will discover that such arcs require 
special handling in our implementation. 

Since a nonsingular F,; must exist, it is well known that the columns of Fi 
represent a rooted spanning tree defined over the nodes of Vg, and that F,, may be 
placed into upper triangular form by row and column permutations. If we choose to 
represent Fi; rather than Fj’ in our implementation, we are interested in performing 


two fundamental operations: 


1. Solving linear systems of the type 


a bi 
and 
23 Fy = ы 


where z; and 2; аге апКпоуп апа 5, and bz are rationals (not necessarily 


integers), and 
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2. Placing Fi, in upper triangular form, where fi; results from a column ex- 
change, a row exchange, a column and row addition or a column and row 


deletion performed on Ё}. 


The literature on (PNE) demonstrates that with a proper choice of data structures, 
operation (1) may be performed very efficiently when b, and b; are integer, and that 
(2) may be done efficiently when the class of updates is limited to column exchanges. 
We will extend these existing approaches to deal with (1) and (2) in the (PNSC) 


setting. 


C. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 
To examine the actions required by Generate.Column(LC), suppose we 
want to generate column LC of the principal part of the tableau (7.3) and place 
= 


the results in the vector z? = (2,,22,23,24)’ which is partitioned to conform to 


(7.3). Rewriting (7.3) in a more convenient manner, we have: 


Б. E 22. 
(1) Б! Ра А) -Fn Fa ~AR En Fa) + Fa -Fū Fi [ÂR (En 7 En Fi! Fill 7 Fi! fis 
0) f -Fn [Án] - Ка |-ЕД! Аа Ау) -Fy -Àn En Fi’) - Fn [Ai (Eis 7. Ei Fi! Fio) 8 Fo 
-Fn -Fp Foi- ÀR En Fn) + Ea] -Fna [ER fas Us < En Fa Ps) ~ Fa Po 
(11) А -Án ЕР АП (ЕЁ = En Fi Fi) 
Gv) | -En AR] - En (- Fn Fai dn) -En -Än En Fi] -En [AR (En - En Fi Fia)| + Ens 
\ -Еп |! А-А Еи КҮ!) + А En [- Fi Fut Ag (Eis 7 Eu Fi! Fia)) Fi! Pis 
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The calculations may then proceed as follows: 
l. Calculate region (iii): 
(a) If column LC is in (7), then: 
z3 = Ae ee 
(b) If column LC is in (77), then: 
z3 = АТВ) 
(c) If column LC is in (777), then: 
z3 = Ат (Ез) = En Fi (Pi) | 


пойсе (һа (Ға) бала ҒІҚҒ3) l may be computed by backpath traversal 


on £F, (e.g., Bradley, Brown and Graves [1977], p. 11). 


2. Calculate region (:): 
(a) If column LC is in (7), then: 
Рүүгү lees 
(b) If column LC is in (77), then: 
Fizi = F223 + e 


(c) If column LC is in (777), then: 


P41 = 1223 + (Бау 
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3. Calculate region (22): 
(a) If column LC is in (7), then: 
za = —F e223 — F121 
(b) If column LC is in (jj), then: 
25 = —Fo223 — £21 
(c) If column LC is in (777), then: 
2› = =Ё›2з. ТЕГЕН 
4. Calculate region (tv): 
(a) If column LC is in (7), then: 
24 = — E2523 — E312 
(b) If column LC is in (77), then: 
24 = — E533 — En 
(с) If column LC is in (777), then: 
= AMAT rS 


Then column LC of (7.4) is available as z? = (z,, 29, 23, 24)’. Note that floating 
point multiplications and/or divisions are necessary only for computation of regions 
(iii) and (iv), and that regions (i), (ii) and (iv) are linear combinations of region 


(iii) and one another. 
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D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSION IN ROW GENERATION 
Now we examine the actions of Generate Row(LR). Suppose we desire to 
generate row LR and place the results of the computation in the vector 2 — (2$, 26, 27) 


which is partitioned to conform to (7.3). Rewriting (7.3) in a convenient form: 


(1) Es - ([- Fi! Fiadi?| En - 1) 55: [- Fi Fudi] Eis» [- ([- Fi! Fo] £i - 1) Fi] Ps 
(ii) | (FaFa Fa- Fan  - {[(Fa Fi Fn - Fadi] Bu > Fan} Fa АЕА = АД) Б. 
+ [- (РАР – РАЈ Е + Р) А) А + Fa 
АЙ, Ат - lái £u fi [А] В + (7 [A7] £u] Po 
(iv) | (En Fn Fa- En) An - ([(En Fi! Еа) А Еа + En] Fi (Е. РР – Еъ) А] Е 
|- (En Fi! fa - En) A] En + Ex} Е Fiy+ Ex | 


Then we may proceed as follows: 


l. Calculate region (7): 


(a) If row LR is in (2), then: 


сы 


5- Sr О) 
(b) If row LR is in (72), then: 


2 (FaFa Fiz = Рз) вА! 


(c) If row LR is in (122), then: 


(d) If row LR is in (tv), then: 
5-(ҒаҒ Fi - En)nÁQ 
2. Calculate region (77): 
(a) If row LR is in (1), then: 
«Ғал-25Ец- бін 
(b) If row LR is in (22), then: 
26 Fy, = —ž5 E11 — (En )LR 
(c) If row LR is in (iti), then: 
Ж«Ғас--ФЕп 
(d) If row LR is in (tv), then: 
Ze Fy, = —2Z5Ey, — (£21) ir 
3. Calculate region (Jjjj): 
(a) If row LR is in (2), then: 
27 = 25 з + 26 Рз 
(b) If row LR is in (22), then: 
1 = Z5 E13 + 56Ё\з + (ЁЮз)рн 
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(c) If row LR is in (222), then: 
27 = 25Ё\з + 26Ё\з 
(d) If row LR is in (iv), then: 
27 = Zg Eja + Ё\з + (ЕЁ›з)гң 


The row LR of (7.5) is available as 2 — (25,26, 27). Note that floating point 
multiplications and/or divisions are avoided for terms involving Ғі2, F?;Fj3 and F3, 


and that regions (jJ) and (jjj) are linear combinations of region (7). 


E. DATA STRUCTURES 

Several suites of data structures have been proposed for algorithms designed 
to solve (PNE) (Johnson [1966], Srinivasan and Thompson [1973], Glover, Karney, 
Klingman and Napier [1974], Bradley, Brown and Graves [1977}). Our implementa- 
tion is based on the last. 

The purpose of the data structures is to define a graph that represents Fi. 
Such a graph forms a rooted spanning tree defined over the nodes of Ng, which 
we denote as Gr,,. Note that because F}, is nonsignular, Fi, must contain at least 
one singleton column. Since Fj; is maintained in triangulated form, we associate 
with each row IR the column JC in which the element appearing on the diagonal of 
IR occurs. This association is analogous to the GUB row/key variable association 
defined for the GUB algorithm. We represent this association with the array KEY(), 
of type INTEGER and dimensioned from 1 to (m +n), which records for each row 
of F,, the extrinsic index of the corresponding key variable, and for each column of 
Е ће extrinsic index of the corresponding basic factored row. KEY() is undefined 


for other row and column indices. 
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We require knowledge of the row ordering of Fiji to define a triangulation. 
PO() is an INTEGER array, dimensioned from 1 to p, which for each row IR of Fy, 
records the extrinsic index of the row that follows IR in the current triangulation. 
This ordering is referred to as preorder, and the successor of IR is called its preorder 
successor. PO() is undefined for factored rows not currently in Ё}. Note that the 
triangulated column ordering of Fi; may be deduced from PO() and KEY(). 

PO() and KEY() capture the triangulated row and column ordering informa- 
tion for a representation of Fi; but provide no means for interpreting this repre- 
sentation as a rooted spanning tree. The predecessor function p() is a well known 
method for representing trees and thus rooted spanning trees. To define p(1,), we 
associate with each node 2; of Gr,, the row index of the offdiagonal element in the 
k'^ column of Fj;, when the rows are ordered to correspond to a triangulation of 
Fii. Graphically, p(z) may be iterated recursively to trace the backpath from node 
i to the distinguished root node of Gr,,. 

We represent p() by the INTEGER array P(), dimensioned from 1 to p, which 
for each row IR in 7i, records the extrinsic index of the predecessor of IR. It is 
convenient to keep a record of the sign of the diagonal element of each row of F1. 
We use the sign bit of P() to do so. Assume node IRP is the predecessor of node 
IR in the current representation of Fi. If the diagonal element in row IR is a plus 
one, then P(/R) = IRP, while if the diagonal element in row IR is a minus one, 
P(IR) = —I RP. In the graph Gr,,, this may be interpreted as follows: if the actual 
orientation of arc a^ = KEY(IR) is the same as that recorded in Gr,,, which is 
(IR, IRP), then P(IR) » 0). If the actual orientation of arc a* is opposite that 
recorded in Gr,,, then P(JR) < 0. 

The final data structure, D(), is an INTEGER array dimensioned from 1 to p. 
For each node IR of Gr,, (row of F\,;) D(IR) records the depth of node IR, where 
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depth is defined as the number of nodes encountered on the backpath between 
IR and the root node. D() allows updates to be performed on the data structures 
representing GF, in a single pass through the data structures, and allows us to avoid 


unnecessary operations in the solution of linear systems to be discussed shortly. 


F. SOLVING LINEAR SYSTEMS 


It is clear from the discussion of Generate Row(LR) and Generate. Col- 


umn(LC) that the solution of the linear systems, 


ШӘ T 
21 Fi = b, 


and 


PFuz-b 


where zı and z, are vectors of unknowns and b, and b, are vectors of rationals, are 
crucial steps. The practical value of our implementation depends to a large extent 
on the efficiency with which these systems may be solved, and thus we will consider 
this problem in some detail. 

First, we examine the data structures used to represent the terms. Three data 
structures are used to represent b, and bz, the right-hand side terms. WORK() is a 
DOUBLE PRECISION array, dimensioned from 1 to p, which is used to hold the 
real values of b, and 5; (we sequence operations so that at any given moment we are 
interested in either b, or b; but not both, so that the same array may be used for 
both). An INTEGER array WORKMK(), also dimensioned from 1 to p, is used as 


a nonzero mask for WORK(). The convention is: 


WORKMK(K) 20 if WORK(K)- 0.0 
WORKMK(K) 21 IF WORK(K) # 0.0 


Finally, the array WORKNZ(), of type INTEGER and also dimensioned from 1 
to p, records the extrinsic column (in bı) or row (in b) index of each nonzero in 
WORK(). The use of WORKMK() and WORKNZ() allows the efficient manage- 
ment of nonzeros in 5, and 5,. 

The solution vectors z; and z2 are represented by three analogous arrays. 
ROWCOL() (for rowcolumn()) is of type DOUBLE PRECISION and is dimensioned 
from 1 to (m+n). It is used to store the representation of a row and column of the 
principal part of the tableau, and thus portions of it are used to store the solutions 
of Equation (7.1). Associated with ROWCOL() is ROWCOLMK(), an INTEGER 
array of conformable dimension, used as a nonzero mask for ROWCOL/(), and the 
INTEGER array ROWCOLNZ(), also of conformable dimension, used to record the 
indices of nonzeros in ROWCOL(). 

In contrast to the usual situation in Simplex-based approaches to linear pro- 
gramming in which each successive right-hand side of the problem may be computed 
by means of a simple update to the previous right-hand side, in this setting such 
is not the case. Thus, we are not able to derive the solution to the current form 
of Equation (7.1) by simply updating the previous solution. Instead, we must solve 
each system from scratch. 

Fach right-hand side is itself the result of a sequence of preliminary computa- 
tions. As each new nonzero element is generated by these computations, its value 


is placed in WORK(), WORKMK() is updated and the index of WORK() in which 


the nonzero appears in placed in the next available position in WORKNZX(). For ex- 
ample, suppose the first nonzero generated for WORK() is the value 6.5 in position 


38. Then the arrays appear as follows: 


WORK(38) = 6.5 
WORKMK(38) - 1.0 
WORKNZ(1) = 38.0 


А counter for the number of nonzeros in WORK() is maintained (an INTEGER 
variable called NNZW), so at the conclusion of the preliminary computations we 
may iterate through the nonzeros on the right-hand side in the order in which they 


were generated by: 


WORK(WORKNZ(K)), K 2 1,..., NNZW 


Since our representation of Fi, is in upper triangular form, the obvious ap- 


proach for solving 


Fuz = b. (1.6) 


is by back substitution. All nonzero elements in F,, are either plus ones or minus 
ones, so only assignments, additions and subtractions are necessary. Assuming that 
the current dimension of Fi, is k by k, we proceed by considering each row IR in 
turn, IR = k,...,1. We assign the value in WORK(IR) to ROWCOL(KEY(IR)) 
and update by an addition the value in WORK() in the row corresponding to the 
predecessor of row IR, (i.e., WORK(P(IR))). 

This approach requires knowledge of the ordering of the rows of Fi, in the 


order of last to first in triangulated form, which is precisely the reverse of our data 
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Figure 7.1: A Triangulated F: 


structure PO(). To support this approach in our implementation, we could define 
and maintain an additional data structure EO(), recording the *endorder" of Р}. К 
we choose not to include an additional data structure, we could instead invert PO() 
in situ whenever this reverse ordering 1s required, inverting it again when PO() is 
next required. With either approach, it is convenient to make PO() and EO(), if 
included, circular lists and to add a distinguished "artificial node", say IRA, which 
is then used to locate both the first node (row) in preorder and the first node in 
endorder. To illustrate, Figure 7.1 displays a triangulated form of Fij with row and 
column labels. 

Figure 7.2 presents the corresponding graph Gr,,. By assuming the existence 
of an artificial node IRA, we in effect change our paradigm of Gr, to that shown in 
Figure 7.3. 

We may interpret the backsolve method as a labeling procedure on Gr,,. In 
this context, each nonzero in the right-hand side is interpreted as a supply or demand 
at the corresponding node in Gr,,. The problem of solving Equation (7.6) is then 


interpreted as one of finding a set of feasible flows defined on the arcs of Gr,,. 
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Figure 7.2: A Basis Graph Gr,, 


The backsolve approach to solving Equation (7.1) has the advantage of sim- 
plicity, and when viewed as a labeling algorithm has the intuitive appeal of “visiting” 
each node only once. The disadvantages are the need for either additional storage 
(EO()) or additional computation (two inversions of PO()). 

The right-hand sides of these problems are invariably very sparse; the density 
is typically 0.3% or less. We are thus motivated to explore alternate approaches to 
solving Equation (7.6). The solution approach just outlined requires the “visiting” of 
every node in Gr,,. Although WORKMK() allows us to perform a simple INTEGER 
comparison to determine if the current node has nonzero supply or demand, we have 
strong empirical evidence that suggests that other approaches are more efficient. It 


is important to keep in mind that this problem is very deeply nested within the 


algorithm structurc, and must be solved tens of tliousands of times for a typical 
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Figure 7.3: An Alternate Graph Paradigm бк, 


linear programming problem. Small changes in efficiency here have large effects on 
the overall algorithm efficiency. 

Since the right-hand sides are sparse, an alternative 15 to consider solving Equa- 
tion (7.6) iteratively as a sequence of subproblems, each consisting of a right-hand 
side containing a single nonzero element and a current partial solution. Iterating 
through the nonzeros, we introduce the supply or demand at the corresponding node 
in Gr,,, and then adjust the flows of all arcs appearing on the backpath of that node 
as necessary. 

Analyzing the worst-case performance of these two approaches, we see that the 
first approach is linear in the number of nodes in Gr,, (rows in F1), while the second 
approach is quadratic in the number of nonzeros in the right-hand side. For the 
densities typically encountered in the problems we have studied, the crossover point 


at which the performance of the linear algorithm overtakes that of the quadratic 15 
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in the range of 400 to 1000 rows in F,,;. However, in the testing we have performed, 
the worst-case analysis is pessimistic, and in practice the performance of the second 
approach always dominates that of the first. 


The second linear system of interest is 


zi Fy = Ы : (7.7) 


We interpret the problem as one of being given flows on the arcs of Gr,, and being 
asked to find the necessary supplies and demands at the nodes. Analogies to each of 
the approaches to solving Equation (7.6) may be found for solving Equation (7.7), 
and our empirical evidence strongly suggests that an approach analogous to the 
second method for solving Equation (7.6) is preferable for solving Equation (7.7). 
Having made the design decisions for solving Equation (7.6) and Equation 
(7.7), it is worthwhile to review our paradigm for Gr,,. Our preferred solution 
techniques do not require a partial ordering of all rows of Fii. Rather, we require 
partial ordering information only among each set of coupled rows, or, graphically, 
among each disjoint component of Gr,,. Rather than introducing an artificial node, 
we may treat each disjoint component as an independent entity. It is then conve- 
nient to treat singleton arcs (either dynamic or static singletons) as self-loops, and 
construct the predecessor function so that it forms a ring within each component. 
Our graph paradigm of Figure 7.1 then becomes as shown in Figure 7.4. The data 


structure representation of Figure 7.4 is then: 


Node: 1 2 3 4 5 6 T 8 
Predecessor: P) 1 -1 1 -4 -4 -5 4 .8 
Preorder Successor: PO() 2 3 1 5 6 Т 4 8 
Depth: D) 0 1 1 0 1 2 1 0 
Key Variable: KEY() 100 101 102 103 104 105 106 107 
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Figure 7.4: The Implementation Paradigm for Gr,, 


G. FACTORED KERNEL UPDATE 

The update of a rooted spanning tree representation of a pure network tri- 
angulated basis which results from a column exchange is well understood and has 
been treated thoroughly in the literature. Bradley, Brown and Graves [1977] give 
an excellent account of such an update in an algorithmic calling very similar to the 
one here. Hence, we will not repeat the details. In this setting, however, a column 
exchange is but one of four possible updates required. We must also deal with row 
exchanges, row and column additions and row and column deletions. Our general 
approach will be to reduce these tliree cases to the column exchange case through 
a sequence of operations which may be thought of as pre- and post-processing. The 
following three cases, along with the column exchange case, comprise the actions 
required for Update_Factored_Kernel. 

The first case we consider is the row and column deletion case, in which the 
dimension of F}; decreases. It is convenient to limit row/column deletions to row/key 


variable pairs. This is because removing a row/key variable pair always preserves 
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nonsingularity of the remaining factored kernel, and because, with our choice of data 
structures, given one member of the pair, it 1s easy to identify the other member 
through the use of the KEY() array. 

With our paradigm for Gr,,, the key variable for row IR corresponds to the arc 
which connects node IR to its predecessor. The removal of IR and KEY (IR) creates 
a new disjoint component within Gr,, consisting of all nodes and the associated key 
variables on whose backpaths IR and KEY(IR) appear. If IR is a leaf (a node which 
is incident to a single arc), then the new disjoint component is null. Within the new 
disjoint component, D() changes for all nodes, but P() and PO() change only for the 
first node in preorder, say IRFIRST, and the last node in preorder, say IRLAST. 
For the first node in preorder, its predecessor becomes itself (with the sign bit used 


to indicate arc orientation), forming a self-loop. The update is thus: 


P (IRFIRST) = IRFIRST 


Its depth is updated to be zero, and the magnitude of this depth change is recorded 
for future use. For every other node in the new disjoint component, D() is reduced by 
the magnitude of the IRFIRST depth change. Finally, PO() of IRLAST is changed 
to the first node in preorder, restoring the local ring structure of PO(). The update 


1S: 


PO (IRLAST) = IRFIRST 


These updates can be accomplished in one pass through the data structures of the 
new disjoint component. The structure of Gr,, is restored by changing PO() of the 
predecessor of IR to the index of the node which was the successor of IRLAST prior 


to the update. 
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The second update case is a row and column addition. Suppose IRADD is the 
row (node) to be added to Fy; (Gr,,). We wish to treat this update as a column 
exchange. To do so, we first incorporate IRADD into the structure and associate 
with it an imaginary (“logical”) arc which forms a self-loop. We then perform a 
column exchange with the new column, say JCADD, entering Gr,, and the logical 
arc leaving Gr,,. 

The difficulty in incorporating IRADD into Gr,, lies with the dynamic single- 
ton arcs. Gr,, currently may have many singleton arcs which appear in our paradigm 
as self-loops. If any of these dynamic singletons are actually incident to IRADD in 
С, incorporating IRADD into Gr,, requires changing their representation. 

To illustrate this point, assume IRADD has the node label “9” and we wish 
to incorporate it into the graph shown in Figure 7.4. Assume also that arcs 103 and 
107 are dynamic singletons which are actually incident to node 9; that is, arc 103 
= (9,4) and arc 107 = (9,8). We first initialize the data structures for node 9 as 


follows: 


р(9) = 0 
Р(9) = 9 
РО(9) = 9, 


which has the effect of placing node 9 at depth 0 as a disjoint component consisting 
of a single node with a (logical) self-loop, as shown in Figure 7.5. 

We then change the representation of arcs 103 and 107 from dynamic singletons to 
doubletons. To do this, we change P(4) from -4 to -9 and P(8) from -8 to -9. We 
then assert a partial ordering among those (formerly) disjoint components which 


have been merged by the change in the representation of the dynamic singletons, 
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Figure 7.5: Initialization of Node 9 Prior to Incorporation into бн 


and retriangulate. This 1s accomplished by increasing the depth of each node in the 
component by one, changing PO(T) from 4 to 8, and PO(8) from 8 to 9. This update 
requires a single pass through the data structures of each affected component. After 
incorporation, Gr,, appears as shown in Figure 7.6. We may now complete the 
update by performing a column exchange. 

To identify the dynamic singletons which are affected by such an update, we 
access the extrinsic problem data structures for row (node) IRADD. For each column 
with a nonzero in row IRADD, we test whether or not that column is currently key 
for some row in Fj. If so, that column is currently represented in Gr,, as a dynamic 
singleton. 

The final update case is the row exchange case. Suppose we want to replace 
row IROUT with row IRIN. We treat this in two stages. The first stage is a row 


and column deletion case, in which the node IROUT and the arc KEY(IROUT) are 


131 


removed from Gr,,. To complete the update, we seek a column (arc) which may 
be added to Gr,, along with row (node) IRIN. We require an arc which is either 
a static singleton, incident to node IRIN, or a doubleton, incident to node IRIN 
and a node (other than IROUT) which is currently in Gp,,. We first consider arc 
KEY(IROUT). If is satisfies the second condition (obviously, it cannot satisfy the 
first condition), we designate it as the arc to be added. If not, we search among the 
variables (arcs) in region (iii) of the tableau for such an arc. We know one must 
exist, for otherwise B is singular. We select the first such arc found as the arc to be 
added. We then perform the row and column addition case. 

Recall that the action Is. Factored Kernel. Singular() is required when we 
have identified an arc (column) for removal from Gr,, and we are considering a 
candidate arc to replace it. We want to determine if this exchange preserves the 
nonsingularity of Fj;. The fact that Gr,, is a rooted spanning tree provides the 
necessary structure to support a simple test of nonsingularity. Assume JCOUT is 
the (known) arc to be removed from Gr,, and JCIN is the candidate replacement. 
In our paradigm, Ff}; is nonsingular if and only if each disjoint component of Gr,, 
is connected and contains exactly one cycle, which must be a self loop occurring at 
the component’s root (the root is the distinguished node in the component whose 
depth is zero). The removal of JCOUT creates a new disjoint component which has 
no root self loop. If the addition of JCIN fails to correct this, the proposed exchange 
is singular. 

The specific test is this: identify the nodes incident to JCIN using the original 
problem data structures. Note that there may be zero, one or two such nodes. If 
none of these nodes are currently included in the structure of Gr,,(Fi1), then the 
proposed exchange is singular. For each such node which is currently included in 


the structure of Gr,,, traverse the backpath of that node by recursively iterating 
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the predecessor array P() until either the root node of that component is found, or 
until the “join” is located (the “join” is that node with largest depth which appears 
on the backpath of both nodes incident to JCIN. The join exists only if JCIN is 
incident to two nodes, both nodes are currently included in the structure of Gr,,, 
and the two nodes are in the same disjoint component.). If JCIN is encountered 
during this backpath traversal, the proposed exchange is nonsingular. Otherwise, it 
Is singular. 

The action Find.Column.to.Remove(IOUT) is particularly simple, since, 
as mentioned previously, we always remove row/key variable pairs. Thus, if IR is 
the node to be removed, then KEY(IR) is the corresponding column to be removed. 

Finally, the action Find.Column.to.Add(IIN) requires searching among 
the variables in region (222) of the tableau. For each candidate arc, we invoke 
Is_Factored_Kernel_Singular(). If the response is “FALSE”, we have found the 


column to be added. Otherwise, we continue the search. 
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VIII. FACTORIZATION OF GENERALIZED 
NETWORK ROWS 


A. INTRODUCTION 


The problem of interest remains: 


min м 
(СМС геи cb F pbyn 
Ly Sr, атоо 
—[y 0; -Inbyn 
where F, E, —I, w, b and r are as before. We now require that each column of 
F have at most two nonzero elements, which may be of arbitrary sign. We may 
associate a generalized network with F by defining a node corresponding to each 


row i of F and an arc F? corresponding to each column j of F. The arcs are defined 


as: 


(7,0) if Fi 0 И О е E 
(0:0) it oO tor all as a 


| | (in ПЕЛҒЫ Ес жапа d s» 
109 0 
We let A= (Fl,..., F^), N 2 (1,...,p) and define the graph G 2 (A/,.A). Arcs 
of the form (2,0) are singleton arcs, sometimes called root arcs. Arcs of the form 
(0,0) are null. 
The most widely studied specialization of (GNSC) is obtained when the F- 


type constraints are equalities and the E-type constraints are vacuous, resulting 


in: 
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пип: шу 
(GNE) 
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Dantzig [1963] provides the seminal treatment of (GNE), identifying the im- 
portant structure that leads directly to efficient primal Simplex-based algorithms 
for its solution. Implementations have been reported by Glover, Klingman, Hultz, 
Karney and Elam [1972, 1973, 1977, 1978, 1979] and Brown and McBride |1984). 

(GNE) may be viewed as a generalization of (PNE), and the cost of such 
generalization is a weakening of the properties which so strongly characterize (PNE) 
and which lead to the efficient and elegant implementations. In (GNE), F is not 
totally unimodular, and thus an implementation may not be restricted to integer 
arithmetic. It is not possible to triangulate every primal Simplex basis extracted 
from F by row and column permutations. The subgraph generated from the columns 
of a primal Simplex basis no longer form a rooted spanning tree defined on the nodes 
of G. Apparently, much has been given up in the generalization. 

In fact, significant structure remains in (GNE). A well known result, due to 
Dantzig [1963], shows that any primal Simplex basis extracted from F can be put in 
the form of Figure (8.1) by row and column permutations. Each square submatrix 
component B* is either upper triangular or nearly upper triangular with only one 
element below the diagonal. Notice that if each B* is strictly upper triangular, the 
structure is analogous to that found in the (PNE) primal Simplex basis. 

Interpreting the structure of B* as a subgraph of G, we find that when B* 
is upper triangular, the subgraph may be viewed as a disjoint component with a 


single self-loop at the root node, exactly as in the (PNE) case. When B* is not 
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upper triangular, the subgraph still forms a disjoint component with a single loop, 
but that loop is no longer a self-loop. For example, Figure (8.2) shows a nearly 


triangular B* with row and column labels. 
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Figure 8.2: A Sample Nearly-Triangulated (GNE) Simplex Basis Com- 
ponent 


The resulting subgraph is shown in Figure (8.3). Such a subgraph is commonly 
called a “one-tree”, (an apparent oxymoron). 

This structure allows the extension of the algorithm and data structures de- 
veloped for (PNE), leading to efficient implementations to solve (GNE) (e.g., Brown 
and McBride |1984)). 

(GNSC) has received less attention in the literature than (PNSC). Hultz 
and Klingman [1976] develop a primal Simplex-based approach to an equality- 


constrained formulation of (GNSC), and report on an implementation which allows 
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a single E-type constraint (Hultz and Klingman (1978)). McBride [1985] develops 
an algorithm for solving a generalization of (GNSC) in which complicating columns 
as well as complicating constraints are permitted, and reports on an implementation 


of that algorithm. 


B. THE FACTORED TABLEAU 

The algebraic development of the prima] row basis B, the nonbasic rows D 
and the factored tableau is exactly as shown in Chapter 7 and is not repeated here. 
Note Һа Д}, the factored kernel, may now be placed in the form shown in Figure 
(8.1) by row and column permutations. The corresponding graph, Gr,,, consists of 
one or more disjoint components, each of which contains exactly one loop. The loop 


may be either a self-loop, as in the case of the disjoint components of (PNSC), or 


— 
W 
= | 


it may be a “root loop” (a loop consisting of two or more nodes, all of which are at 


depth zero in the component), as shown in Figure (8.3). 


C. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN COLUMN GENERATION 
The sequencing of computations for column generation in (GNSC) is exactly 


the same as that for (PNSC). 


D. SEQUENCING COMPUTATIONS TO EXPLOIT COMMON SUB- 
EXPRESSIONS IN ROW GENERATION 
The sequencing of computations for row generation 1n (GNSC) is exactly the 


same as that for (PNSC). 


E. DATA STRUCTURES 

Several proposed data structures for algorithms tailored to solve (GNE) have 
been offered (Glover, Klingman and Stutz [1973, 1974], Elam, Glover and Klingman 
[1979], Brown and McBride [1984]). Our implementation is based on the last, which 
extends to (GNE) the data structures developed in Bradley, Brown and Graves 
[1977] to solve (PNE). 

We have seen the similarity between the structure of the disjoint components 
that arise in the graph corresponding to Fi, in (GNSC) and the structure of the 
components in (PNSC). To develop a representation of these components, we wish to 
extend the data structures developed for (PNSC) in a natural way. Knowledge of the 
row ordering of Fj; is again maintained in an INTEGER array PO(), dimensioned 
from 1 to p. As before, it is convenient to make PO() a circular list defined on 
each disjoint component. The unique correspondence between each row of Fy, and 


a distinguished column is maintained in the INTEGER array KEY(), dimensioned 
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from 1 to (m -- n). The conventions for KEY() are exactly as in (PNSC). The depth 
array, D(), is defined exactly as before. However, as suggested by our diagram 
in Figure (8.3), we adopt the convention of placing all nodes appearing on a loop 
at depth zero, which is consistent with our definition of *root loop". Finally, a 
representation of the predecessor function is needed, and we define P() as before. 
For any row IR in the factored kernel, P(IR) » 0 indicates that the diagonal element 


in the near-triangulation is the first non-zero factored element in its column. 


F. SOLVING LINEAR SYSTEMS 


Just as in (PNSC), the crucial steps in the generation of rows and columns of 


the principal part of the tableau are solving the systems 


NM = pi (8.1) 


and 


Fiz-b (8.2) 


where z; and z3 are vectors of unknowns and b, and b; are vectors of rationals. Our 
approach for solving these systems closely parallels that developed for (PNSC). 
The data structures used to represent 21,22,06, апа 8; аге ехас у the same 


as in the (PNSC) implementation. To represent 6, and 62, we use: 


WORK() 
WORKMK() 


WORKNZ(), 
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and for z, and 22: 


ROWCOL() 
ROWCOLMK() 


ROWCOLNZ(). 


The right-hand side sparsity-exploiting approach for solving (8.1) and (8.2) 
developed in Chapter 7 may still be used, but two complications arise in the (GNSC) 
setting. First, the nonzero elements of Fj, are no longer restricted to be plus and 
minus ones, and may assume arbitrary values. Thus, we are not able to restrict all 
arithmetic operations to additions and subtractions. 

Interpreting (8.1) and (8.2) as problems of finding balancing supply and de- 
mands (8.1) or feasible flows (8.2), we see that in (GNSC) we must allow for gains 
and losses in the flows over the arcs of Gr,,. Because our formulation allows two arbi- 
trary nonzeros in each column of F (rather than one arbitrary nonzero and one plus 
one, for example), backward and forward substitution schemes require both multi- 
plications and divisions. We may eliminate the need for divisions by pre-computing 
both ratios of the nonzero elements in each arc. That is, if a column has nonzero 


elements а and b in rows of F, we pre-compute the ratios ? and 2. We then substi- 


tute the pre-computed value for the division whenever it occurs. Storing a pair of 

DOUBLE PRECISION real numbers for every column of (GNSC) is wasteful, and 

instead we define a single pointer for each column which points to the location of 

the first of two DOUBLE PRECISION real values. The first is the value ¢ and the 
b 


second is 7. (Note that for singleton arcs, the values a and + are stored instead.) 


We compute and store only the unique ratio pairs for a given problem instance, and 
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the number of such unique pairs is usually quite small: typically, a problem with 
10,000 to 20,000 columns has only 30 to 40 unique ratio pairs. 

The second difficulty is that the disjoint components of Fj; may be nearly tri- 
angular rather than triangular, and thus backward and forward substitution cannot 
be applied directly. However, Dantzig [1963] shows that a variation of backward 
and forward substitution may be used to solve (8.1) and (8.2). This method solves 
the triangulated part of a disjoint component exactly as is done in backward and 
forward substitution. When a root loop is encountered, the method requires two 
calculations for each row or column of Fi; corresponding to a node or arc on the 
loop. The second calculation involves a term Dantzig [1963] called the *loop factor", 
which is common value for every node (or arc) on the loop. Our implementation 
uses this modified forward and backward solution technique for solving (8.1) and 
(8.2), and we store unique values of loop factors in à DOUBLE PRECISION array 
called GNROOT(). 


G. FACTORED KERNEL UPDATE 

We first consider the actions required by Update. Factored Kernel. 

Brown and McBride [1984] give an excellent treatment of the update required 
for a column exchange within a generalized network basis in an algorithmic setting 
very similar to the one here, and thus we do not repeat that discussion. As in 
the case of (PNSC), however, three additional classes of updates may occur in this 
setting: row exchanges, column and row additions and column and row deletions. 
We again treat these three cases by reducing them to the column exchange case 
through sequences of pre- and post-processing operations. 

We again limit row and column deletions to row/key variable pairs. When 


row IR and its associated key variable KEY(IR) are removed from fF}, a disjoint 
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component is created which has no loop. If IR and KEY(IR) do not appear on a root 
loop of Gr,,, a new disjoint component is created, just as in the (PNSC) algorithm. 
If IR and KEY(IR) do appear on a root loop of Gr,, the resulting retriangulated 
component contains a self-loop rather than a root loop. The update of the subgraph 
data structures is very similar to those presented in Chapter 7. 

Row and column additions require the incorporation of the incoming node into 
the structure of Gg, prior to the column exchange, just as in (PNSC). The technique 
for (GNSC) is exactly the same as in (PNSC). Once the new node has been added 
to GR, a column exchange operation is performed, with a logical arc (which forms 
a self-loop on the new node being added to Gr,,) being replaced by the new arc 
(column) being added. 

Finally, row exchanges are handled as a two-stage update, exactly as in (PNSC). 

The task of determining the singularity of F1; (Is -.Factored Kernel. Singu- 
lar) is more challenging in (GNSC) than in either of the other two algorithms due 
to the arbitrary nonzero elements. In each of the previous two algorithms, we have 
been able to determine with certainty the singularity of Fi; indirectly. In the (GUB) 
factorization, Ғіі = Aj, a signed identity matrix. Since A), is orthogonal (i.e., 
Ai, - Ai, = 1), Fir is perfectly conditioned (see e.g., Golub and Van Loan [1983] for 
a discussion of matrix conditioning). In (PNSC), F1; may be triangulated, placing it 
in a form with plus and minus ones on the diagonal. Thus, the determinant of Fj; is 
either plus or minus one. We may therefore determine its singularity by examining 
the structure of Cg, rather than considering Fi; directly. In either case, as long 
as the accumulated round-off error in the current representation of the problem 
Is such that we can distinguish plus or minus one from zero, we may rely on the 
structure of our representation of Cg, to deduce the singularity of Fj,. Thus, we 


are able to discern singularity logically and need not resort to analytic methods. 
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In (GNSC), however, we may not rely exclusively on the structure of Gg, to reach 
conclusions about the singularity of Fj,. It is not difficult to invent examples of 
factored kernels corresponding to instances of (GNSC) which are ill-conditioned. 
We expect ill-conditioning of Fj, to lead to serious numerical difficulties, which 
we wish to avoid if possible. Since we generally have freedom in determining the 
structure of Fi, and since our fundamental rank-one update of F3; is the column 
exchange, we use a heuristic based on column exchanges to anticipate conditioning 
problems. Recall in the column exchange update, a column has been selected for 
removal from Fi; and a candidate column has been identified as its replacement. 
We use the backpath traversal method described in Chapter 7 to determine if the 
proposed exchange maintains the required structure of Gp,,. If it does not, we can 
conclude that the exchange renders Fi, singular and we reject the candidate arc. 

As we traverse the backpath, we perform calculations which, if the exchange 
does in fact preserve the required structure of Gr,,, will ultimately compute the 
determinant of the disjoint component in which the replacement arc will appear if the 
proposed exchange is performed. If the absolute value of the computed determinant 
is too small (a decision controlled by a implementation parameter), we conclude that 
the proposed update produces a submatrix of Fy}; (the submatrix corresponding to 
the new disjoint component) that may be ill-conditioned, and thus Fi; may itself 
be ill-conditioned. We therefore reject the proposed candidate arc. 

This approach is attractive because it is computationally inexpensive and, 
based on our empirical evidence, it works well in practice. It is of course merely 
a heuristic, since it is easily shown (e.g., Golub and Van Loan [1983)] that the 
determinant can be a poor indicator of matrix condition. Further, we do not actually 
compute the determinant of F,,, but rather only the determinant of a submatrix of 


Fa. 
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The remaining update actions, Find.Column.to. Add(ICI) and Find_Col- 


umn_to_Remove(ICI) are treated in a manner analogous to that in (РМ5С). 
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IX. COMPUTATIONAL RESULTS 


A. INTRODUCTION 

The algorithms developed here have been implemented and extensively tested 
within the framework of a commercial-quality optimization system: the X-System 
of Brown and Graves [1975]. This system employs the Graves mutual primal-dual 
algorithm in a variety of large scale optimization applications, including linear, non- 
linear, mixed integer and decomposed models. Although we report computational 
results only for linear models, our factorizations are seamlessly integrated into the 


X-System and support all system features. 


B. TEST PROBLEMS 

The benchmark test suite is drawn from a wide variety of actual applications. 
Table (9.1) provides a short synopsis of each model, quoting from the abstract 
where a reference in the open literature is available. In some cases, several different 
instances of models are reported. We selected these models because they provide 
a representative sample of size, structure and taxonomy of contemporary real-life 
applications of linear programming. For those models that are mixed integer, we 


report solution statistics for the initial linear programming relaxation. 


TABLE 9.1: Description of Test Suite Models 


e GTE “The seven Telephone Operating Companies within GTE have adopted 
an integrated business system called Capital Program Management System 
(CPMS) to guide their 3 billion dollar per year capital planning. The system 


includes a large scale mixed integer programming optimization system that 


145 


optimizes the critical economic tradeoffs between maximizing the long-term 
budget value of the firm's equity and satisfying shorter-term financial con- 
straints, resource limitations and service objectives. Investment opportunities 
for the next 5 years are modeled as 0-1 variables with alternative implemen- 
tations for each. The objective is to maximize the net present value of the 
capital portfolio. There are financial constraints on capital, internally gen- 
erated funds, net income to common, and limits on resources such as labor 
hours, lines installed, etc. There are also constraints that enforce logical re- 
lationships among opportunities (such as, if choose A then must choose B).” 


See Bradley [1986]. 


INVEST “Capital allocation and project selection for a large multi-national 
firm is modeled as a two-stage multi-year nonlinear capital budgeting problem 
with over 40,000 integer variables. Innovative modeling yields subproblems 
easy to solve, and optimality is achieved with a single iteration of the nonlinear 
master problem.” See Harrison, Bradley and Brown [1989]. The instance 


reported here is a linear program subproblem of the two-stage model. 


TANKER “A crude oil tanker scheduling problem faced by a major oil com- 
pany is solved using an elastic set partitioning model. The model takes into 
account all fleet cost components, including the opportunity cost of ship time, 
port and canal charges, demurrage and bunker fuel. The model determines 
optimal speeds for the ships and the best routing of ballast (empty) legs, as 
well as which cargoes to load on controlled ships and which to spot charter. All 
feasible schedules are generated, the cost of each is accurately determined and 


the best set of schedules is selected.” See Brown, Graves and Ronen [1987]. 
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e GAS "Natural gas utilities supply about one fourth of the energy needs of 
the United States. From wellhead to consumer, operations are governed by an 
astounding diversity of purchase, transport and storage contract agreements 
which enable complex physical distribution systems to meet future demands 
no more predictable than next year’s weather. Gas is a highly detailed op- 
timization model which utilities use to plan operations and to justify such 


plans to regulatory agencies is developed.” See Avery, Brown, Rosenkranz 


and Wood [1989]. 


e ODSM “А commonly occurring problem in distribution system design is the 
optimal location of intermediate distribution facilities between plants and cus- 
tomers. À multicommodity capacitated single-period version of this problem 
is formulated as a mixed integer linear program. A solution technique based 
on Benders Decomposition is developed. ... An essentially optimal solution 
is found and proven with a surprisingly small number of Benders cuts.” See 
Geoffrion and Graves [1974]. The instances reported here are decomposition 


master problems. 


e TAM “The annual decision on how much of the Air Force procurement budget 
should be spent on the many different aircraft and how much should be spent 
on the many different munitions is of great interest to many people. How 
the Air Force staff develops information to support the decision has changed 
over the years. Currently, a linear program is being used by the Air Force 
Center for Studies and Analysis and is being tested by the Munitions Division 
of the Plans and Operations Directorate (AF/XOXFM) for munitions tradeoff 
analysis. The LP uses existing data and estimates on (1) aircraft and munition 


effectiveness, (2) target value, (3) attrition, (4) aircraft and munition costs, 
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and (5) existing inventories of aircraft and munitions. Other factors considered 


are weather and length of the conflict." See Might [1987] and Jackson [1989]. 


PHOENIX *The U. S. Army operates a fleet of over 7,000 helicopters to 
perform combat and combat support tasks. Although newer, more technically 
advanced helicopters have been and are being procured, the majority of the 
fleet is still composed of helicopters that were built during the late 1960s 
for use in South Vietnam. These older airframes are rapidly reaching the 
end of their useful lives and must be (i) replaced by newer, more advanced 
designs, (ii) gutted and refit or replaced by a combination of (i) and (ii). Army 
force planners recognize that this problem can only be solved by a long-term 
program and the commitment of billions of dollars. Phoenix is a software 
system employing mixed integer linear programming to help Army aviation 
staff and officers develop long- range helicopter modernization plans.” See 


Clemence, Teufert, Brown and Wood [1988]. 


AMMO4H “А four-commodity transshipment model for delivery over time of 
military products from production and storage locations to overseas locations 
to support theater operations is developed. The model covers five physical 
echelons, including production plants, storage depots, ports of embarkation, 
ports of debarkation and geographic field locations. Road, rail, sea and air 
transportation are modeled, and product demands are time- phased. Capaci- 
tation occurs primarily on sea and air links, and as throughput capacities on 
transfer points, requiring replication of some echelons. The objective of the 


model is to minimize deviation from on-time deliveries." See Staniec [1984]. 


GK A weekly multi-plant production/inventory/transshipment linear pro- 


gram from a consumer products industry is developed. The model is meant 
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to guide weekly processing and packaging decisions. Production consists of 
two stages: basic products are produced and then packaged into different- 
sized containers to yield finished products. Processing lines typically produce 
a subset of the basic products and have limited capacity with overtime charges 
for weekend shifts. Packaging lines for finished products are similar. In-house 
inventory capacity is limited although outside storage is available at additional 


cost. Inter-plant shipments are made by rail or truck. See Wood [1989]. 


C. METHODOLOGY 

We wish to evaluate implementations of our three algorithms on the basis 
of computation time and computer memory requirements. Since each algorithm is 
simplex-based, the formal theory of algorithmic complexity provides no basis for 
preferring one to another, since in the worst case none enjoys a measure of running 
time that 1s polynomial in the size of the problem specification (see e.g., Garey 
and Johnson, [1979] for a discussion of algorithmic complexity and Klee and Minty 
[1972] for an analysis of the Simplex method). Thus, we are led to consider “typical” 
performance by gathering empirical evidence on the performance of the algorithms 
solving “typical” problems. 

We prefer an implementation that is both fast and requires little computer 
memory. Most researchers who have reported on implementations of related al- 
gorithms have been concerned primarily with execution speed, and certainly it is 
important. However, we have seen in our algorithmic setting that once high speed 
storage has been allocated for the problem representation and for the program code, 
all remaining memory is available to store the representation of the explicit trans- 
formation kernel and the factored kernel. If the solution trajectory is such that their 


combined size never exceeds available memory, we succeed in solving the problem. 
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If not, we fail. When success or failure depends on the total storage requirements 
of the inverse representation, we may be willing to sacrifice execution speed in ex- 
change for an economical representation of the inverse. This is a classic theme in 
computer science and software engineering, and we believe its importance in this 
context has been largely overlooked. 

We will be comparing the performance of four separate implementations. The 
first is an unadorned version of the X-System, which implements the Graves mutual 
primal-dual method as presented in Chapter 2. There is no underlying factorization. 
We refer to this implementation as ^X". The second implementation is the GUB 
factorization presented in Chapter 6, and referred to as "GUB". The third is the 
pure network factorization of Chapter 7, referred to as “PNET”, and the last is the 
generalized network factorization, called “GNET”, of Chapter 8. 

The ideal approach for this computational study, would be to develop four 
equivalent formulations of each model, each customized for its particular implemen- 
tation with the goal of inducing a large factored row set of the appropriate type. 
This approach is a consistent theme in the literature dealing with specialized al- 
gorithms and one that we strongly endorse. Alternate formulations of a model are 
often available, and it seems sensible to choose one that as strongly as possible 
exploits the strengths of the solver. 

However, all of the models used here are “off-the-shelf” in the sense that they 
were developed at various times by various modelers, and we cannot afford to develop 
alternate formulations. Thus, our approach is to preserve a single, unfactored rep- 
resentation of each model, and attempt to identify favorable row structures through 
the use of heuristics. Our procedure is based on the work of Brown and Thomen 
[1980], Brown and Wright [1983] and Brown, McBride and Wood [1985]. The heuris- 


tics are greedy and myopic in the sense that they initially consider the entire row 
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set of the problem, and discard one row at a time without backtracking until the 
remaining set satisfies the desired row factorization. This may easily result in the 
confounding or destruction of structure intended or perceived by the modeler. While 
our approach yields interesting and useful observations about the implementations, 
it is in some ways a poor substitute for the customized model method. 

Table (9.2) tabulates the important structural information concerning the 
model instances we will be solving. The column headings may be interpreted as 
follows: m is the total number of structural constraints (note that this use is differ- 
ent from that in the problem specifications of Chapters 6, 7, and 8), n the number 
of variables, pgyg the number of GUB rows identified by the heuristic, ppy the 
number of pure network rows, pgn the number of generalized network rows and NZ 
the number of nonzeros in the technological coefficient matrix of the model. For 
example, consider the first model in Table (9.2), GTE. The structural constraints 
contain 57,563 nonzeros, and the model consists of 6,624 variables. When viewed 
as an unfactored mutual primal-dual model, it consists of 960 explicit constraints 
and 0 factored constraints. When viewed as a GUB factorization, it consists of 909 
factored (GUB) constraints and 960 —909 = 51 explicit constraints. Similarly, when 
viewed as a pure network factorization (PNET), it contains 909 factored rows and 51 
explicit rows. Finally, when viewed as a generalized network factorization (GNET), 


it consists of 922 factored rows and 38 explicit rows. 


D. COMPUTATIONAL RESULTS 

We solved each of these problem instances on an IBM 3033/AP under the 
VM/CMS operating system using VS FORTRAN 1.4.1. A virtual machine size 
of six megabytes was used, simply because it is the largest size normally available 


to us. It is the nature of a time-shared system that measurements of processing 


times are somewhat imprecise due to system load factors and accounting techniques 
for system processing overhead. We have attempted to mitigate these effects by 
performing our experiments during periods of low system usage. 

Table (9.3) displays the solution times for each of the test problem instances by 
each of the four implementations. Àn **" indicates that the solver failed to solve the 
problem instance. Such a failure occurred because the storage requirements for the 
explicit transformation kernel representation exceeded the memory available and the 
solver terminated gracefully. Solution times represent only the CPU time required 
to solve the problem, and exclude the initial problem input, the time required by 
the factorization heuristic to identify the factored row structure and the final output 
to record the solution. All figures are 1n CPU seconds. 

The formulations of three of the test problems were strongly influenced by the 
design of the target solver: the TANKER model possesses a strong GUB structure 
since it contains a set of schedule selection constraints for each ship (1.e., from a set 
of candidate schedules, select at most one). The AMMOA4H model is a multicom- 
modity capacitated transshipment problem and is thus best suited to a pure network 
factorization. The PHOENIX10 model design was shaped by the generalized net- 
work factorization paradigm. The nature of the factored row structures shown in 
Table (9.2) supports this assertion. In the TANKER model, the factored pure net- 
work rows are exactly the same as the the factored GUB rows, and the heuristic 
constructs a generalized network factored row set by identifying one additional row 
to be paired with each GUB row. The AMMO4H model may be viewed as a GUB 
factorization with a relatively modest GUB set consisting of the joint capacitation 
constraints, or as a pure network factorization with a relatively large pure network 
(PN) factored row set. In the PHOENIX10 model, the dominant structure is clearly 


the generalized network row structure. 
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As we expect, the factorization is most successful when the model is wed to the 
solver. Although we are surprised to find the performance of (PNET) competitive 
with (GUB) on the TANKER model, (GUB) clearly dominates the X and GNET 
solvers. Similarly, (GNET) dominates on the PHOENIX10 model and (PNET) on 
the AMMO4H model. We would be disappointed if the results were otherwise. 

We see evidence in Table (9.3) to suggest that the approach of using heuristics 
to automatically identify factored structure has its pitfalls. In a number of problem 
instances, although our heuristics identify significantly larger factored sets as we 
progress from the base system to (GNET), we see little improvement in computation 
times (INVEST, ODSM1, TAM8). In fact, we see in the TANKER model that 
the temptation to confound the modeler’s intended GUB structure by specifying a 
generalized network factorization leads to disastrous consequences, even though this 
tactic doubles the size of the factored row set. Interestingly enough, when we apply 
the (GNET) solver to the (GUB) factored row set, we are able to solve the problem 
in 16.5 CPU seconds. This suggests that the “quality” of a row factorization is not 
completely specified by size alone. 

We are encouraged by the observation that the transition from the basic system 
to (GUB) to (PNET) seldom degrades solution times, even when doing so yields 
little gain in the number of additional factored rows. This seems to contradict 
popular folklore, which suggests that computation times worsen as the transition 
from unadorned Simplex to (GUB) to (PNET) is made unless the transition is 
accompanied by a substantial increase in the size of the factored portion of the 
model. In fact, computational testing reported by others is frequently limited to 
models in which the number of explicit rows is in the range of one to twenty (see 
e.g., Chen and Saigal [1977], Glover, Karney, Klingman and Russell [1978], Glover 


and Klingman [1981]). Our results are all the more remarkable given the the lack of 
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guidance from the modeler for the *intended" row factorization. We note, however, 
that results are mixed for the transition from (PNET) to (GNET), and it seems 
clear that the applicability of the (GNET) factorization is not as general as that of 
(PNET). 

Finally, we observe that in several models it is factorization that separates 
success from failure in solving the problem with a given allocation of computer 
resources. This fact alone may be reason to consider this approach in practice. 

Our second interest with respect to computation 1s in the memory requirements 
of our algorithms. Our design strategy allocates memory to the data structures 
which represent Fj, so that we may successfully represent the largest dimension of 
factored kernel that may possibly arise. Limited memory remains to store the repre- 
sentation of Az, the explicit transformation kernel. During the solution process, if 
we encounter a representation of Aj] which requires more memory than is available, 
failure occurs. We measure the size of the computer representation of Az} and the 
amount of available computer storage in terms of the elements of Aj} that can be 
stored. The number of bytes per element varies according to the size of the problem 
(this has to do with the FORTRAN data types INTEGER*2 and INTEGER*4), but 
is generally 28 bytes per element. Table (9.4) lists for each problem instance/solver 
pair the maximum size of the Aj} representation encountered during the solution, 
measured in units of number of elements. An asterisk (*) indicates that the number 
shown equals the maximum number of units of storage that were available, and thus 
failure occurred. 

We see that generally the representation of the maximum size of the explicit 
transformation kernel decreases as the generality of the factorization increases. Re- 


calling the definition of the explicit transformation kernel: 
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this trend is as we would expect. As the generality of the factorization increases, 
we expect the size of the factored component to increase and the size of the ex- 
plicit component to decrease. Each potentially binding explicit row which may be 
converted to a factored row by adopting a more general factorization reduces the 
dimension of the resulting representation of Aj. Also, the density of the term 
— Fy, Fyy’ Fi2 generally increases as the complexity of the structure of Fy; increases. 
Assuming the dimension of Fj; is k by k, the number of nonzeros in Fi! in the GUB 
factorization is k. In the pure network factorization, the number of nonzeros in F1! 


2 


may be as large as + 


>» and in the generalized network factorization, the number of 


nonzeros in Fj;! may be as large as n?. We note that (GNET) again provides several 
exceptions to the general trend in Table (9.4). 

It is the dynamic nature of our factorization algorithms which marks this 
work as a departure from previous research. Table 9.5 illustrates the significance 
of this point. The first column lists the number of constraints which are binding 
at optimality, and the second column expresses this as a percentage of the total 
number of constraints in the problem instance. The results shown here are typical 
of real-world large-scale models. It is usually the case that many constraints are 
not binding at optimality, and there are computational advantages to be gained by 
exploiting this fact. | 

Columns 3, 4 and 5 of Table 9.5 list for (GUB), (PNET) and (GNET) respec- 
tively the number of explicit constraints that are binding at optimality. Since in 
each implementation, binding factored constraints are handled more efficiently than 


binding explicit constraints, we see that the computational success of our dynamic 
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factorization algorithms is due to the fact that even in large model instances, we 
are able to limit our attention to a relatively small number of explicit constraints, 
usually on the order of a few hundred or less. While this is well beyond the size of 
previously reported implementations, our results show that it is quite manageable. 

It is useful to establish a criteria for comparing and contrasting the perfor- 
mance of the algorithms which accounts for both execution time and storage re- 
quirements. Although more sophisticated models could undoubtedly be developed, 
we offer a simple model which we feel captures the essential features we wish to 


consider. Define 


FSM з! 


where s is the total computer storage (measured in megabytes) required to solve 
a problem instance (including program code, original problem data, tableau rep- 
resentation, factored kernel representation and explicit transformation kernel rep- 
resentation) and ¢ is the execution time (measured in CPU seconds). f(s,t) is 
monotonically increasing in s and t, and in a crude fashion it captures the essential 
features of the way in which computer resources are often marketed commercially. 
We will use f(s,t) as a measure of performance. Table (9.6) displays f(s,t) for each 


иж? 


of the problem instance/solver combinations. indicates that the solver failed to 
solve the particular instance. 

We observe that the trend as the transitions are made from base system to 
(PNET) is a decline in f(s,t). It is apparent that (PNET) is a versatile imple- 


mentation, performing extremely well on models with highly favorable structure 


(AMMO4H, KG4, GASPNA, GASPNC) and comparing favorably with (GUB) and 
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(GNET) on every other model. The performance of (GNET) is generally quite corn- 
parable to that of (PNET). It appears to be about 15-20% slower than (PNET) 
when it is used to solve a pure network factorization (GASPNC, AMMO4H). Our 
row elimination heuristics are usually effective in identifying a favorable factored 
structure, but the computational results on the TANKER model clearly illustrate 
the limitations of this approach. (GNET) apparently requires a more sophisticated 
and careful user than do (GUB) or (PNET), and in this sense it is perhaps a more 
specialized algorithm. The evidence shown here indicates that (PNET) is a strong 
candidate for use as a general implementation, and need not be viewed as a highly 
specialized implementation suitable only for rare model instances. For the models 
studied here, we have progressed well beyond the stage of solving instances with 


only a handful of explicit rows. 
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TABLE 9.2: Summary of Problem Suite Dimensions 


m 
GTE 960 
INVEST 1,338 
TANKER 93 
GAS PN A 6,848 
GAS ENIC 3,194 
GAS PNE 1,184 
GAS GN A 6,848 
GAS GN C 3,974 
GK 2 3,819 
GK 3 5,128 
GK 4 1,636 
ODSMI 3,023 
ODSM2 594 
TAMI 91 
TAM2 180 
TAM3 209 
ТАМ4 211 
ТАМ5 438 
ТАМ8 420 
TAMI? 629 
PHOENIXIO 1,618 
PHOENIX30 4,305 
AMMO 4H 13,963 


n 


6,624 
11,989 
7,598 
27,884 
15,362 
5,102 
27,884 
15,362 
17,844 
27,493 
37,139 
11,568 
ООДА 
389 
1,204 
2.883 
i07 
10,969 
6,104 
17,793 
6,884 
4,291 
83,497 


PGUB 


909 
941 
33 
4,345 
2.658 
434 
4,484 
2.664 
1,265 
2.295 
2.498 
523 
490 
22 

42 

63 

59 
102 
118 
177 
206 
293 
1,071 
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PPN 


909 
1,101 
33 
5.934 
3,418 
877 
5,142 
3,084 
2,578 
3,867 
5,156 
540 
490 
28 

54 

81 

T 
132 
154 
231 
220 
303 
12,892 


PGN 


917 
1,168 
66 
5,976 
3,490 
882 
5,976 
3,420 
2,585 
3,876 
5,167 
558 
490 
34 

66 

99 

98 
162 
196 
295 
1,153 
3,605 
12,892 


М2 


57,563 
36,829 
30,890 
36,702 
19,701 

7,355 
36,702 
19,701 
34,809 
54,280 
73,766 
21,522 
42.821 

1,212 

6,869 
21,356 

6,954 
93,964 
49,376 


164,947 


13,818 
16,441 


166,713 


TABLE 9.3: Solution Times in CPU Seconds 


X СОВ РМЕТ СМЕТ 
СТЕ 48.4 33.0 37.2 38.3 
INVEST 21.6 17.1 il ried 
TANKER 141.8 20.5 17.7 43.8 
GASPNA à КІШ 2271 
GAS PN C 96.9 93.7 35.8 31.3 
GAS PN E i 146.8 17.4 20.1 
GAS GN A ^ б * 5216 
GAS GN C 2229 50.9 37.6 33.7 
GK 2 27.1 226 2/153 19.4 
(519 62.4 58.4 50.4 45.1 
GK 4 : 380.9 182.2 186.6 
ОП5МІ Son 231 9.3 8.8 
ODSM?2 34.1 32.5 32.8 30.4 
ТАМ1 0. 0.7 0.7 0.7 
TAM2 2.5 ki 22 2 
TAM3 17.8 14.1 14.3 13.5 
ТАМ4 1.2 IE 1.5 ШЕТІ 
ТАМ5 317.1 3000 307.8 269.6 
ТАМ8 93.3 81.8 83.8 77.4 
ТАМ12 660.9 660.6 624.1 480.9 
PHOENIX10 68.6 43.7 30.4 10.7 
PHOENIX30 1 1 * 495.1 
AMMO 4H МА 11,2413 2530 288.4 


e NÀ indicates problem instance not run. 


t This problem was run on an IBM 3081K under the MVS operating system using 
VS FORTRAN 1.4.1 in a 32 megabyte virtual machine using an advanced starting 
solution. The solution time shown is adjusted to account for the approximate dif- 


ference in computing speed between the IBM 3081K and the IBM 3033/AP. 
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СТЕ 
INVEST 


TANKER 

GAS PN A 
GAS PN C 
GAS PNE 
GAS GANTA 
GAS GN C 


СК 2 
GK 3 
GK 4 
ODSMI 
ODSM2 
ТАМІ 
ТАМ2 
ТАМЗ 
ТАМ4 
ТАМ5 
ТАМ8 
ТАМ12 


РНОЕМІХ10 
PHOENIX30 
AMMO 4H 


X 


11,571 
8,380 
2.256 

*132,995 

21,912 
*174,482 
%132,995 

17,966 
8,280 

13,991 

“110,334 
1,329 
418 
[Т 
3,042 

10,430 
1,489 

29,661 

19,545 

37,716 

Т ІН 

*153,781 
NA 


GUB 


22] 

682 

415 
*130,883 
14,765 
111,821 
*130,883 
8,017 
2,719 

D o 
19,646 
861 

3 

1.231 
2,078 
8,226 
1,202 
23,190 
14,886 
28,935 
38,705 
*151,746 
122,085 
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PNET 


214 
654 

S 
35,395 
260 
4,808 

* 127,780 
1,258 
221 

388 
6,529 
259 

3 

841 
1,444 
5,305 
867 
15,689 
8,351 
16,097 
21,879 

* 148,348 
23 


TABLE 9.4: Number of Elements in Explicit Пааво Е Kernel Rep- 
resentation of Optimality or at Failure 


GNET 


196 
521 

81 
65,558 
364 
4,219 
65,081 
608 
184 
310 
9,155 
421 


1,165 
1,610 
6,595 
1,239 
18,737 
11,182 
20,486 
1,458 
11,523 
23 


TABLE 9.5: Summary of the Number of Binding Explicit Constraints at 
Optimality 


Total Percent GUB PNET GNET 


Binding 
GTE 592 RE 19 18 18 
INVEST 158 56.7 199 194 162 
TANKER 50 60.2 ШІ 30 9 
GASPNA 3,051 44.6 t 292 349 
CAS PN C 2-3 61.6 849 66 9] 
GASPNE 897 ОТОТ 88 81 
GAS GN A 3,045 44.5 1 t 352 
GAS GN C ЖЕ 58.1 863 402 88 
GK 2 1,950 ОРО 93 90 
GK 3 2,912 50.5 1,260 140 141 
СК 4 4,119 53.9 216 363 356 
ODSMI 297 9,8 53 49 47 
ODSM2 448 15:4. ЛО 215 0 
TAMI 46 50.5 44 32 35 
ТАМ? 81 48.2 ТТ 5] 61 
ТАМЗ 148 БЕ ШЕШІ 97 107 
ТАМ4 121 Bi LUG 59 81 
ТАМ5 25 SL OD 170 186 
TAMS 216 ы zs 140 174 
TAMI?2 413 65. RS 208 252 
PHOENIX10 1,085 ӨТІ 1084-7122 TT 
PHOENIX30 3,477 80.8 t : 109 


AMMO 4H 2,889 20.7 2,882 7 7 
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TABLE 9.6: f(s,t) in Megabyte - 


X GUB 
GTE 84 4T 
INVEST 25 16 
TANKER 112 15 
GASPN A * t 
GAS PN C 73 58 
GAS PNE * 500 
GASGN A $ * 
GAS CNC 62 45 
CK? 27 19 
GK 3 90 74 
GK 4 + Sees 
ODSMI 4 6 
ODSM2 49 47 
TAMI 0.1 0.1 
ТАМ? 1 1 
ТАМЗ 14 10 
ТАМ4 0.1 0.1 
ТАМ5 150 662 
ТАМ8 133 106 
ТАМ19 29707 2.015 
PHOENIXIO 169 65 
PHOENIX30 t г 
АММОАН МА 9,093 


РМЕТ 


54 
19 
13 
926 
21 
9 

* 
30 
18 
60 
309 
5 
48 
0.1 


0.1 
614 
95 
1,882 
32 


* 


933 


seconds 


GNET 


64 
20 
42 
А 
31 

14 
Wt 
34 

21 
66 
387 
7 

53 
0.1 

1 

19 
0.1 
622 
110 
1,632 


693 
1194 


X. CONCLUSIONS 


We have presented three dynamic row factorization algorithms for solving 
large-scale linear programs. Although each may be used to solve any LP instance, 
each is designed to exploit a particular model row structure: generalized upper 
bounds (GUB), pure network rows (PNET) or generalized network rows (GNET). 

Previous research by others generally suggests that specialized algorithms such 
as those presented here are useful only when the factored structure completely dom- 
inates the structure of the model instance. There are reports of algorithms for 
solving problems having a single unfactored (explicit) constraint (Hultz and Kling- 
man [1978], Klingman and Russell [1978]). When implementations are reported, 
problem suites are limited to instances having a very small number of explicit con- 
straints, typically in the range from one to twenty (Chen and Saigal [1977], Glover, 
Karney, Klingman and Russell [1978], Glover and Klingman [1981]). The consensus 
seems to be that such algorithms are appropriately viewed as specialized algorithms, 
useful only for solving very special problem instances. 

Our experience strongly refutes this view. We find the performances of our im- 
plementations of the dynamic factorization algorithms are competitive with that of 
a commercial-quality optimization system on every model instance we have tested. 
This is particularly remarkable for two reasons. First, our test suite consists of 
models developed by skilled modelers specifically to exploit the capabilities and 
characteristics of the solver with which our implementations are competing. Sec- 


ond, we must select the row factorizations without the benefit of guidance from the 
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modeler, relying instead on useful but imperfect heuristics. Despite these computa- 
tional handicaps, our tests show our implementations to be at least as efficient as a 
well-respected commercial-quality optimization system. 

Our development has stressed the similarity between the algorithms and the 
natural extension which leads from one to the next. This is in contrast to the devel- 
opment that has been reported for similar, non-dynamic algorithms (e.g., Dantzig 
and Van Slyke [1967], Klingman and Russell [1978] and Hultz and Klingman [1978]) 
in which the specifics of the individual algorithm obscure the generality of the ap- 
proach. The conceptual difference between our algorithms is seen to be largely 
isolated to the structure of a single algebraic entity, the factored kernel. By abstract- 
ing the structure of the factored kernel and concentrating on the general algorithm 
design, we demonstrate the versatility and flexibility of this approach. 

We are gratified to find that the modularity suggested by the algorithmic 
development can be realized in an implementation design. We succeed in developing 
a software suite which displays a “single-system image”. The modularity of the 
algorithm allows the definition of an “abstract data type” (see, e.g., Aho, Hopcroft 
and Ullman {1974]) which isolates the data structures and update procedures for the 
factored kernel from the rest of the implementation. Each factorization 1s seamlessly 
integrated within the system design, presenting a single design image. 

The early 1980’s produced a great deal of research in the area of automatic 
identification of special structure in LP models (see, e.g., Gunawardane and Schrage 
[1977], Glover [1980], Schrage {1981], Brown, McBride and Wood [1985] and Bixby 
and Fourer [1986]). We have incorporated the most useful of these ideas into our 
implementation, and we have what we believe to be the first implementation which 


supports the automatic identification of factored row sets. This capability may 


164 


be used to identify new factored structure or to validate or augment a modeler- 
provided recommendation. Our computational experience indicates that while this 
approach is not as promising as perhaps first envisioned, it is nonetheless a valuable 
tool. When faced with the choice of either solving an unfactored model instance 
or automatically identifying a factored structure and then using the corresponding 
solver, our results show that the latter is nearly always to be preferred. Our results 
seem to suggest, however, that in addition to quantity of factored rows, the issue 
of quality of factored rows exerts influence on the performance of the factorization 
algorithms. While not well understood, it is clear that the myopic approach of 
our heuristics is no substitute for the modeler’s guidance in identifying factored 
structure. 

Several areas suggest themselves for further research. Certainly additional fac- 
tored structures can be examined. For example, one approach for treating factored 
column structures (“complicating columns”) is to allow the partitioning of rows into 
the categories “factored” and “explicit” to vary as the algorithm progresses. That 
is, allow factored row set membership to be determined with respect to the column 
structure of the currently nonbasic variables rather than with respect to the column 
structure of all problem variables. While conceptually simple, such a generalization 
seems to present significant algorithmic challenges. 

General algorithms are sometimes useful in specialized contexts. For example, 
processing networks (Koene [1982]) are network models which allow proportional 
flow restrictions on the arcs entering or leaving some nodes. One formulation of such 
a model results in a pure or generalized network structure with a set of complicating 
columns. Chen and Engquist [1986] propose a primal partitioning algorithm for 


solving processing network problems. An alternate formulation yields a pure or 


generalized structure with ОШОЛ rows, and we note that this is precisely the 
structure we seek for our network factorizations. 

The multicommodity capacitated transshipment problem (MCTP) has been 
the subject of much research over the years, and a number of specialized algorithms 
(see, e.g., Assad [1978] or Kennington [1978]) have been proposed to solve it. Adopt- 
ing a general perspective, MCTP may be viewed either as a GUB model or as a pure 
network model with side constraints, and either view might be preferred depending 
upon the dimensions of the particular instance under consideration. Our compu- 
tational experience indicates that the pure network factorization algorithm offers a 
powerful technique for solving MCTP. As an experiment, we customized our (PN) 
implementation to exploit the special structure of the side constraints in MCTP. It 
is interesting to note that in our scientific computing environment, we observed no 
difference in solution times between the customized version of (PN) and the original 
implementation. 

Finally, all the approaches we have considered assume the prior existence of a 
specific structure in the factored rows which in turn determines the structure of the 
factored kernel. An extension of this general approach is to relax the requirement for 
strict conformance to a specific structure. Instead, we might allow the factored row 
structure to be “nearly” homogeneous. For example, we may allow a small number 
of complicating columns to disturb what is otherwise a factored pure network row 
structure. We then expect the structure of the factored kernel to be dominated by 
that induced by the predominant row structure, with only occasional complications 
due to the exceptional row structure. We allow for this exceptional structure in 
the factored kernel by identifying it “on-the-fly” as the algorithm progresses, and 
treating it in an appropriate manner. This approach may be thought of as a hy- 


brid between the factored method developed here and dynamic basis triangulation 
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methods (see, e.g., Hellerman and Rarick [1971] and [1972], Saunders [1976] and 
McBride |1980)). 

Dynamic extrinsic factorization is subsumed if we activate functions in the 
update analogous to the secondary exchanges now employed. Essentially all that 
has to be done is ensure that successive factored components retain their stipulated 
special structure. In our estimation, this will only be justified in cases where the 
model structure is amenable, and quite likely will require some model-specific fea- 
tures to perform well on difficult models. We have limited our experimentation to 


those static extrinsic cases which are believed to be most useful. 
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